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ABSTRACT 


This report describes the development of an analytic procedure for the 
calculation of nonequilibrium boundary layer flows over surfaces of arbitrary 
catalycities. A n existing equilibrium boundary layer integral matrix code was 
extended to include nonequilibrivm chemistry while retaining all of the general 
boundary condition features built into the original code. For particular appli- 
cation to the pitch-plane of shuttle type vehicles an approximate procedure was 
developed to estimate the nonequilibrium and nonisentropic state at the edge of 
the boundary layer. 

The nonequi librium code (BLIMF/KINET) was used to calculate catalycities 
of typical shuttle thermal protection materials which were exposed to an arc 
jet environment. Sufficient data was available to predict the simultaneous 
catalytic efficiencies for atomic oxygen and atomic nitrogen recombination. 

These catalycities and appropriate nonequi librium edge conditions were used to 
predict the boundary layer behavior on the pitch-plane of a typical shuttle 
vehicle over a trajectory range which included Loth laminar and turbulent flows. 
These calculations show a small reduction in heat transfer when compared to 
catalytic surfaces; however, substantial reductions are demonstrated for 
noncata.lytic surfaces. These calculations also demonstrate the significance of 
including entropy layer effects. 
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SYMBOLS 


parameter used in the solution of the mixing length equation 
(defined by Equation (126) ) 

radial offset of the generator axis at angle of attack. Figure 6-4 

caloric enthalpy deviation associated with equilibrium dissociation 

parameter used in the solution of the mixing length equation (de- 
fined by Equation (127)) 

coefficients in polynomial body surface equation 
binormal direction unit vector, body oriented coordinates 


constant introduced in the a constraint (Equation (44) ) 

constant introduced in the approximation for multicomponent thermal 
diffusion coefficients embodied in Equation (8) . Tentatively 
established by correlation of data to be -0.5 


product of density and viscosity normalized by their reference 
values (defined by Equation (52)) 

tangent shock surface metric. Equation 178 
frozen specific heat of the gas mixture 

property of the gas mixture whi'h reduces to when diffusion 

coefficients are assumed equal for all species (defined by 
Equation (8) ) 

specific heat of species i 
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D 




V 



ERROR 
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coefficients defined in finite-difference representation of 
streamwise derivatives (defined in Equations (107) and (108) for 
two- and three-point difference relations, respectively) 


a reference binary diffusion coefficient introduced by the approx- 
imation for binary diffusion coefficients embodied in Equation (7) 

square of ratio of tangent to normal shock surface metrics, 

(C/N ) 2 

multicomponent thermal diffusion coefficient for species i 
multicomponent diffusion coefficient for species i and j 
diffusion coefficient for all species when all V 1 are equal 

binary diffusion coefficient for species i and j 
unit direction cosines, jth principal axis 

errors for the various equations during Newton-Raphson iteration 
(driven toward zero in the iteration) 

stream function (defined by Equation (45) ) 

diffusion factor for species i introduced by the approximation 
for binary diffusion coefficients embodied in Equation (7) 

ratio of effective polytropic exponents in equation of state, 

Y/(Y-1) 

Gibbs free energy fox j-th specie 

static enthalpy of the gas (def inod by Equation (5) ) 


Vi 


SYMBOLS 

(.continued) 


static temperature 


velocity component parallel to body surface 


local .dimensionless tangential velocity component 


shear velocity, defined in Equation (33) 


free stream velocity 


velocity component normal to body surface 


local dimensionless normal velov ty component 


mole fraction of species i 


dimensionless axial measure with origin at nose intercept of 
body generator axis 


dimensionless axial trace of shock wave with origin at generator 
axis intercept 


truncated series obtained in Taylor series expansion of 


/ f'p dn (defined by Equation (112)) 


distance from surface into the boundary layer, measured normal 
to the surface 


dimensionless normal measure above body surface 


dimensionless y-coordinate defined by Equation (33) ) 
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SYMBOLS 

(continued) 

total number of elements; also nixing length constnat 


mass; fraction of molecular species i 


total mass fraction of element (or base gas) k contained in sur- 
face material (e,g, , charl removed by combustion, sublimation, or 
vaporization 


total mass fraction of element (or base gas> k contained in gas 
which enters boundary layer without phase change at the surface 
(e.g,, pyrolysis gases) 


total mass fraction of element (or base gas) k irrespective of 
molecular configuration (defined by Equation (11) ) 

partial pressure equilibrium constant for m-th chemical reaction 

mixing length (defined by Equation (39) ) 

dimensionless mixing length (defined by Equation (57) ) 

parameter used in mixing length formulation (defined by Equation 
(123)) 


mass flow rate per unit area 


mass removal rate per unit area of surface material (e.g., char) 
by combustion, sublimation f or vaporization 


dimensionless mass flow associated with the ith streamline 


mass flow rate per unit area of gas which enters boundary layer 
without phase change at the surface (e.g., pyrolysis gases) 

mass removal rate per unit area of Jl** 1 component surface material 
(e.g., silica) in the condensed phase (e.g., by melting with 
subsequent liquid runoff or by spallation) 


molecular weight of the gas mixture 
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molecular weight of species i 
Mach number 

effective shock Mach number, Equation 162 - 

number of nodal points across the boundary layer selected for the 
purpose of the numerical solution procedure 

normal direction unit vector, bod” oriented coordinates 

normal shock surface metric. Equation 178 

dummy variable representing f', H^, or 

dimensionless pressure, ith region 

dimensionless total pressure 

pressure, also a parameter used in the mixing length formulation 
(defined by Equation (120)) 

partial pressure of species i 

frozen Prandtl number of the gas mixture (defined by Equation (71) ) 

turbulent Prandtl number (defined by Equation (55) ) 

diffusional heat flux per unit area away from the surface 

heat conduction per unit area into the surface material 

one -dimens iona 1 radiant heat flux (toward the surface), that is, 
the net rate per unit area at which radiant energy is transferred 
across a plane in the boundary layer parallel to the surface 
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metric coefficient for streamline spreading (equal to lcc-tl r idiut. 
in the boundary layer in a meridian plane for axisymmetrxc flow) 


surface value of r 


dimensionless radial displacement measured from generator axis 


universal gas constant 

Reynolds number; subscripted with the length scale if other than s 
effective nose radius for Newtonian flow 
reaction rate for m-th reaction (Equation 91) 
nose radius of curvature 

distance along body from stagnation point or leading edge 

initial streamline arc measure. Equation 199 

reference system Schmidt number (defined by Equation (74)) 

turbulent Schmidt number (defined by Equation (54)) 

radial offset of forward stagnation point from generator axis. 
Figure 6-4 

mass fraction of i-th specie 

parameter de' iued to simplify problems with transverse curvature; 
see Equation ll)) 

tangential direction unit vector, body oriented coordinates 
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static temperature 


velocity component parallel to body surface 


local dimensionless tangential velocity component 


shear velocity, defined in Equation (33) 


free stream velocity 


velocity component normal to body surface 


local dimensionless normal velov ty component 


mole fraction of species i 


dimensionless axial measure with origin at nose intercept of 
body generator axis 


dimensionless axial trace of shock wave with origin at generator 
axis intercept 


truncated series obtained in Taylor series expansion of 


/ f'p dr) (defined by Equation (112)) 


distance from surface into the boundary layer, measured normal 
to the surface 


dimensionless normal measure above body surface 


dimensionless y-coordinate defined by Equation (33) ) 
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constant in the mixing length, differential equation (see Equation 
(.29) ) 


Y principal normal coordinate of body 


Z principal binormal coordinate of body 


Z^ a quantity for species i which is introduced as a result of the 

approximation for binary diffusion coefficients and reduces to K. 
when all diffusion coefficients are assumed equal (defined by 
Equation (8) ) 

A 

Z angle of attack annular radius of shock wave surface. Figure 6-4 


a quantity for element (or base species) k which is introduced as 
a result of he approximation for binary diffusion coefficients 
and reduces to when all diffusion coefficients are assumed 
equal (defined by Equation (12)) 


ZP^, ZP 2 «.*. truncated series obtained in Taylor series expansion of integrals 
S involving nonsimilar terms (defined by Equation (118)) 

v~ angle of attack of body generator axis 


a* flux normalizing parameter (defined by Equation (67) ) 




normalizing parameter used in definition of n (see Equation (43)) 
defined implicitly by use of a constraint such as Equation (46) 

mass fraction of clement ( or base species) k in species i 

streamwise pressure-gradient parameter (defined by Equation (53)) 

effective poly tropic exponent, equation of state 

catalytic efficiency ^ 
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y-diniension normalizing parameter (defined by Equation (56)) i 
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6 . 

x 


dimensionless shock standoff at ith station 


l A l-l 


logarithmic distance between two streaowise positions denoted by 
the subscripts l and t*-l (defined by Equation (109)) 


Af ^ ,Af| corrections for f^ , f^, . . « , during Newton-Raphson iteration 


AlS® 


change in free energy for j-th chemical reaction 


6 * 


6* 

1 


displacement thickness (defined by Equation (36)) 


incompressible or velocity displacement thickness (defined by 
Equation (37) ) 


6n 


distance between two boundary layer nodal points (defined by 
Equation (100) ) 


inverse density jump at shock front j-th station 


A 

n,n 


transformed coordinate in a direction normal to the surface (de- 
fined by Equation (47)). Note: the hat is dropped from n through- 

out most of the report 


angle between a surface normal and a normal to the body center- 
line; also time in discussions of the charring ablation program 


0 


angle measured clockwise from normal to body generator axis. 
Figure 6-3 


thermal conductivity 


dimensionless curvature 


y 


shear viscosity 
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properties of the gas mixture (defined by Equation (8)) 
which reduce to unity, to M, to 1/M, and to An M, respectively, 
for assumed equal diffusion coefficients 


stoichiometric coefficient of j-th reactant in m-th chemical 
reaction (Equation 90) 


stoichiometric coefficient of j-th product in m-th chemical 
reaction (Equation 90) 

kinematic viscosity 

dimensionless streamline distance, origin at shock intercept 

transformed streamwise coordinate (defined by Equation (47)). 
Mote: the hat is dropped from r throughout most of the report 

density 

dimensionless density 

total mass flux per unit area into the boundary layer 


individual species turbulent eddy diffusivity 


average turbulent eddy diffusivity, where it is assumed that all 


P £ r 


P £ r 


turbulent eddy conductivity 


turbulent eddy viscosity 


dimensionless eddy viscosity (defined by Equation (59)) 
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Stefan-Boltzmann constant 


local shear stress 


third body efficiency of i-th specie in m-th chemical reaction 


elemental production rate (Equation 12) 


rate of mass generation of species i per unit volume due to 
chemical reaction 


Subscripts 


refers to baseline trajectory case, a 
Altitude = 164 Kft 


30°, M « 9.61. 

OO 


refers to a catalytic surface 


pertains to boundary- layer edge 


equil 


pertains to surface equilibrium requirement 


pertains to the i** 1 species or to the i^ 1 nodal point in tho 
boundary layer, starting with i = 1 at the surface 


pertains to j species 


pertains to k element (or base species) 


pertains to SL streamwise position 


pertains to m iteration during the Newton-Raphson iteration 
process 


pert? ' o to the n" nodal points , corresponding to the outer 
edge ot the boundary layer solution 
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SYMBOLS 

(concluded) 


generator axis conditions 


pertains to the stagnation point 


pertains to the steady state energy balance requirement 


pertains to wall 


reference condition, usually taken as zero streamline from 
inviscid solution (synonymous with boundary-layer edge in the 
absence of an entropy layer) 


Superscri pts 


refers to conditions at shock wave 


equal to unity for axisymmetric bodies and zero for two-dimensional 
bodies 


signifies that quantity is normalized by a* (e.g., j* = 

/V 

represents partial differentiation with respect to T) or n (usually 
tj unless otherwise noted) 
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SECTION 1 
INTRODUCTION 

One of the basic features of the space shuttle vehicle is the reusable thermal 
protection system which, depending upon the number of flights between refurbishments , 
strongly affects the economics of transporting men and equipment into space. Be- 
cause of the restrictions imposed by current high temperature materials technology, 
equilibrium radiation surface temperatures must be held to relatively low values 
compared to surface temperatures on a typical RV. For the space shuttle vehicle, 
low surface temperatures will be achieved by aerodynamic deceleration at high 
altitudes. However, the low densities associated with high altitudes raises 
questions about the effects of nonequilibrium chemistry and the validity of bound- 
ary layer assumptions, especially in the stagnation region. The size of the pro- 
posed shu' le vehicle and its complex geometry poses questions about the effects 
cf entropy gradients and shock interference heating. In addition the designers 
must be concerned about transition and turbulent flow, gap heating, and dis- 
continuities in stir face chemistry due to the tile-like assembly of the TPS. 

This reprrt describes the development of analytic procedures to calculate 
the effects ! nonequilibrium chemistry in laminar and turbulent boundary layers 
over sux. ces of arbitrary catalycities . In order to effectively address this 
problem, however, it was also necessary to establish a procedure for estimating 
t' •. chemical state of the air at the edge of the boundary layer. This state is 
determined from the inviscid streamline history which therefore also accounts for 
envxopy layers . 

The computational procedure has been used to predict surface heating rates 
01 tl i windward pitch plane of a typical shuttle vehicle and was compared with the 
more normal equilibrium boundary layer predictions. The procedure has also been 
used to examine the effects of chemical reaction rates, surface catalycity and 
stagnation pressure as well as to predict catalycities of candidate shuttle TPS 
materials in simulated reentry environments. 
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SECTION 2 

NONEQUILIBRIUM BOUNDARY LAYER CODE - BLIMP/KINET 

The development of the governing boundary layer equations and computational 
procedure for the simultaneous solution of these equations were presented in 
Reference 1. The procedure described in Reference 1 is applicable to nonsimilar 
multicomponent laminar boundary layers with arbitrary equilibrium or nonequilibrium 
chemistry, unequal concentration and thermal diffusion, radiation absorption and 
emission, second order transverse curvature effects, and a genert set of surface 
boundary conditions which includes an intimate coupling with transient charring- 
ablation energy and mass balances. A turbulent eddy viscosity model is presented 
in Reference 2. The resulting Fortran IV computer code which was developed in 
accordance with the analyses of References 1 and 2 is designated by the code 
name BUMP (Boundary Layer Integral Matrix Procedure) but does not include the 
radiation-emission model and, although permitting selected surface rate-controlled 
reactions, is restricted to equilibrium chemistry in the boundary layer. The 
radiation version is code named RABLE and is described in Reference 3. In the 
following, the extension of BLIMP to include nonequilibrium chemistry in the 
boundary layer is discussed, (Kinetically controlled surface reactions were al- 
ready included in the BLIMP code as an input option,) 

Since the analysis and computational procedure has been presented in 
References 1 and 2 the following will be only a brief reiteration of these 
references with an emphasis on those features which are germaine to homogeneous 
chemistry and surface catalyzed rate controlled reactions. It should be emphasized 
that none of the generalized procedures and capabilities developed in BLIMP have been 
destroyed or altered in the process of extending BLIMP to BLIMP/KINET. BLIMP/KINET 
is currently limited to 7 nonequilibrium species controlled by a maximum of 20 
stoichiometric chemical reactions, exclusive of different third body efficiencies. 
Kinetic rates are calculated in accordance with the development in Reference 4. 
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2.1 MATHEMATICAL MODEL 

2.1.1 Conservation Equations 

Omitting radiation transport, the global conservation of mass, momentum and 
energy, and the specie conservation equations are respectively 
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and the coordinates r, s, and y are shown in Figure 2-1. The flux for multi- 
component diffusion is obtained from the Stefen-Maxwell relation 



( 6 ) 


using the bifurcation approximation of References 5 and 6. In this procedure, 
the binary diffusion coefficient is approximated by the function 



D(T,P) 

F i F j 


(7) 


where D is a reference diffusion coefficient and is a diffusion factor for 
specie 1. Then with a few of the following definitions 
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the Stefen-Maxwell equations can be solved explicitly for the diffusive flux, i.e., 


j = - K + (z _ K ) 
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In addition the diffusive energy flux can be expressed as 
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For some problems (e.g., equilibrium chemistry) the number of differential 
equations to be solved can be substantially reduced if conservation of "elements" 
rather then conservation of species is used. Thus by defining the elemental mass 
fraction by 
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the diffusive flux will be represented by 
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Since elements are conserved, E a.i, « 0 so that this formulation would not be 

j 1 j j 

appropriate for nonequilibrium chemistry since the production terms vanish. How- 
ever by reintroducing the production term into this conservation equation, the 
same formulation can be used for equilibrium or nonequilibrium chemistry by the 
simple expedient of setting the production terms to zero for equilibrium condi- 
tions. In addition, for nonequilibrium chemistry each specie is considered as 
an element and the "elemental*' specie conservation equation serves to alter the 
elemental composition of the gas mixture. The conservation equations to be 
solved are therefore 

global (15) 
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These equations, which include the approximations for unequal thermal and 
multicomponent diffusion coefficients of Reference 6, are parabolic in nature and 
thereiure requiring specifications of the dependent variables, their derivatives, 
or a linear combination thereof along the wall (y = 0) , edge of the boundary 
layer, and at the initial body station. Typical sets of boundary conditions will 
be discussed later in this report. Also necessary in the mathematical formula- 
tion of the problem is the specification of the molecular transport properties, 
equation of state and equilibrium relations for the multicomponent gas, and a 
description of the eddy viscosity, conductivity and diffusivity. The molecular 
transport properties, equation of state, and equilibrium relations are discussed 
in references 2 and 4. 

2.2 TURBnLENT FLOW CONSIDERATIONS 

In tha conservation equations developed above, the concepts of eddy 
viscosity, eddy diffusivity, and eddy conductivity were used to express the 
correlations of fluctuating velocity, species, and enthalpy fields in terms of 
mean field quantities. This is only one of several possible techniques of closing 
the set of equations (assuming satisfactory expressions for the eddy parameters 
are available) , and it does not provide any information regarding the evolution 
of the turbulent correlations as the flow progresses downstream. Admittedly, it 
would be more desirable to describe the turbulent fluctuations in a more complete 
manner such as with an entrainment relation, turbulent kinetic energy relation, 
or a local turbulent constitutive equation (Reference 7) . However, thece tech- 
niques are still in early stages of development even for incompressible single 
component flows, therefore a more proven approach was selected for the present 
analysis. The Boussinesq description of turbulent boundary layers has proved to 
be very useful, particularly for complex reacting flows such as are being de- 
scribed here, and will be used exclusively in the present analysis. 

There is a wide amount of latitude possible even within the eddy viscosity 
framework of turbulence, particularly in applying classical incompressible 
models to compressible flows. The following two subsections describe how the 
turbulence model described in Reference 8 was applied to the present compressible 
flow problem. 
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a. Wall Region 


Following the work of Clauser (Reference 9) the boundary layer is 
divided into a law of the wall region and a wake region. The relatively thin 
wall region of the turbulent boundary layer is characterized by very steep 
gradients in the turbulent transport and mean field properties. Turbulent 
stress varies from zero at the wall to near its maximum value the outer edge 
of the wall region. There is a vast amount of empirical e ice that these 
turbulent stresses and also the mean flow field propertiet be described en- 
tirely in terms of the wall state, wall fluxes, thermodynaii »;d trsuisporc pro- 

perties of the fluid, and the normal coordinate y. Since the streamwise coordin- 
ate does not enter the solution for this region, the problem becomes a one- 
dimensional initial value problem. Eliminating s derivatives from the continuity 
equation and neglecting variations in r due to the thinness of the layer results 
in 


o 

dy 


(19) 


or 


pv ** p v 
w w 


( 20 ) 


where the subscript w refers to the wall value. Thus the wall injection rat , 
^w V w' w *’* c * 1 he a function of s, determines the transverse mass flux through 
the entire wall region. Using the same technique for the momentum equation 
and substituting equation (20) 
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K w w M dy w 

where the wall shear, T w , is also typically a function of s. For flows over an 
impermeable wall with constant properties, this equation reduces to 
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or 

T - T w (23) 

indicating that shear can be considered constant in the wall region. For flows 
w:th injection or ablation, it is seen that shear varies with the mass injection 
rate and local velocity, that is, 


T » T + p v u (24) 

WWW 

This one-dimensional description of turbulence in the wall region will be useful 
in formulating a mixing length model for eddy viscosity as described in the fol- 
lowing paragraphs. It should be made clear however, that only the wall region 
turbulent shear stress is assumed to behave in a one-dimensional fashion. In 
the solution procedure, the complete two-dimensional equations of motion are 
solved over the entire boundary layer. 

A complete investigation of the validity of the mixing length postulate 
for flows with injection has been reported in Reference 10. The analysis used 
in this investigation is an extension of that work ; therefore, the reader should 
refer to Reference 10 for more details. 

Because of the current lack of understanding of turbulent mechanisms, 
"theoretical" predictions of the variation of turbulence near the wall must rely 
on empirical input into relations based on some phenomenological dependence. The 
generality of the ultimate goals of this analysis and the desire to approximate 
the physical situation dictated certain prerequisites for the turbulent transport 



relations. 

These were: 

[ 

a) 

The relationr must indicate a continuous variation of the turbulent 
transport properties from the wall to the fully turbulent region. 

4 

* 

t 

i 

i 

b) 

The relations must be generally applicable to mass, momentum, and 
energy transport. 

y> 

r 

r 

c) 

The relations must be applicable to compressible or incompressible 
flows with real gas propc-ttes. 


d) 

The relations should be suitable for transpired and untranspired 
boundary layers without any, or a minimum, modification cf form. 


I 
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Two basic variations of the eddy viscosity hypothesis have been proposed 
in the past. The first type predicts the variation of turbulent viscosity from 
the wall to the fully turbulent region. The second type of hypothesis involves 
a variation of mixing length from the wall into the fully turbulent portion of the 
boundary layer. Data indicate that surface mass addition strongly affects the 
eddy viscosity profile, and it was found that tht first type of hypothesis could 
not be simply iwodified to predict this variation. On the other hand, success 
of the mi> . ng length .theory in predicting profiles in the fully turbulent por- 
tion of the boundary layer with surface mass addition has been noted, for ex- 
ample, in Reference 11 and 12. It has generally been concluded that the slope 
of the linear relation between mixing length and distance from the wall is in- 
sensitive to surface mass addition . As a consequence of this apparent generality 
of the mixing length approach, it vas adopted for the present studies. 

The basic mixing length postulate can be expressed as 

(ov1 ' ’ p*’ {%)’ - <*„ 37 (25 > 

where the mixing length, t, is a combination of various correlations, but retains 
some relationship to the scale of turbulence. Prandtl proposed that this length 
will, in its simplest form, be related to the distance from a wall, at least in 
the region of development of turbulence. His proposition that 


d* . . _ 

— « constant, K 
dy 


(26) 


has been tested under a variety of conditions and found to be quite adequate in 
the fully turbulent portion of the wall region. 

As the wall is approached however, this simple relation is no longer 
appropriate, and, in fact, it can be shown theoretically that 


iim l - n j 

y-*o w » 

. . dJt . ( 

lxm dy “ 0 \ 

y-»o * J 


( 27 ) 
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This is a consequence of the Reichardt-Elrod criterion (see Reference 10). Thus, 
two criteria are specified, namely, Prandtl's hypothesis which is appropriate in 
the fully turbulent portion of the wall region and the Reichardt-Elrod wall cri- 
terion as expressed by Equation (27) . 

Several means of expressing a relation covering the full range of y and in- 
cluding these limiting criteria have been used by other investigators. It is 
advantageous in considering extensions of mixing length theory to establish some 
physical logic for the selected relation. Unfortunately, the understanding of 
transition from the laminar to the turbulent portions of the layer has not reached 
a state permitting any quantitative specification. Therefore, the selected model 
can be based only on qualitative understanding of the process, dimensional con- 
siderations, and the above limiting criteria. These criteria are satisfied for 
incompressible flows by a simple implicit relation of the form 

~ « (Ky - *) (28) 

which implies that the r\te of increase of the mixing length is proportional to 
the difference between the value postulated by Prandtl (Ky) and its actual value. 
This rate of increase is assumed to be augmented by the local shear and retarded 
by the local viscosity. Using these parameters to nondimensionalize the above 
relation yields 


dSL 

ay “ 


(Ky - £) 



(29) 


where y is the constant of proportionality. The' coefficients K and y were 
& & 

shown in Reference 10 to be invariant for a wide variety of flow conditions at 
values of 0.44 and 11.83, respectively. 

For compressible flows, the physical arguments must be changed somewhat. 
Rather than describing the scale of a turbulent eddy, it seems appropriate to 
describe the mass of the eddy, pZ, with respect to the mass available, f P dy. 
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Thus, by analogy to Equation (29), the rate of increase of the mass of an eddy 
will be taken to be proportional to the difference between the mass available 
between the wall and the point of interest (times an appropriate constant) and 
the mass of the eddy: 



(30) 


Nondimensionalizing as above. 


1 dg£ 

P dy 



(31) 


The constants K and y + are left at their incompressible values of 0.44 and 11.83 
for the time being. The integral-differential character of this mixing length 
equation indicates a difficult solution procedure in the physical coordinate 
plane. However, in the (n,5) coordinates introduced by the Levy-Lees transfor- 
mation, the mixing length equation simplifies somewhat. This will be discussed 
further in Section 2.3. 

For the special case of constant properties and zero injection (constant 
shear) , Equation (31) can be integrated to yield 


l 





(32) 


where 



+ 

y 


yu T 


v 


(33) 
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It can be seen that the Reichardt-Elrod criteria is satisfied at the wall. For 
large y, the expression 



is obtained. This special case result for constant property zero injection flows 
is not used in the general analysis technique presented here. 

b. Wake Region 

The wake region of a turbulent boundary layer is so named because the 
flow in this region tends to have a wake-like character. In particular, the 
outer 80 to 90 percent of the boundary layer combined with the local turbulent 
eddies dominates the mixing processes within the flow, and the viscous effects 
become second order. Gradients in the wake region are typically much smaller 
them those of the wall region. Since the pressure gradient and streamwise deriv- 
ative terms are important in the wake region, the two-dimensional character of 
the tutbulence must be considered in its entirety, as opposed to the approximations 
of the wall region, 

A fortunate feature of the wake portion of the boundary layer is that eddy 
viscosity is nearly constant across this region, at least for equilibriumt in- 
compressible flows. In particular, Clauser (Reference 9) was able to relate the 
eddy viscosity to edge velocity and a length scale 6* 

e„ = 0.018 u.6* (35) 

for a great quantity of experimental data taken in equilibrium flows. 

The quantity 6 * in this relation is the displacement thickness 


6 * 



( 36 ) 


'Equilibrium as used here refers to a particular pressure gradient, (d*/T w ) 
(dP/dx) , which results in self-similar velocity profiles (Reference 9) . 
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in which the densities cancel out £or incompressible flows. For compressible 
flows, this length scale is inappropriate since under some conditions 6* can be 
negative. Defining a velocity defect thickness as 



( 37 , 


; the eddy viscosity in the wake portion o f the flow will be taken as 

e„ ~ 0.018 u.6 * (38) 

H XX 

i 

< a satisfactory technique for choosing the correct £ M expression at any particu- 

*, lar body station is to use the wall region expression 

V 

\ s. ■ f < 39 > 

% 

until e exceeds the wake value. Equation (38), at which point e is held con- 
stant at the wake value for the remainder of the boundary layer thickness. 

c. Boundary Layer Transition 

As can be seen from the form of the conservation equation, both the 
> molecular and turbulent transport terms are considered simultaneously. This is 

• necessary since an accurate description of the turbulent boundary layer requires 

• that the time-averaged fluctuation terms disappear near the wall. Another reason 
for the inclusion of these terms is the description of laminar or transitional 

• flows. From the form of Equation (38), it can be seen that for very small 6* 

| the turbulent stresses will be small compared to the laminar ones. Without any 

l constraints on the equations as stated above, kinematic and eddy viscosities are 

9 . 

| equal at a velocity displacement thickness Reynolds number of 56: 

T 

e u 0.018 u,5* 

H XX 

V V 

■ 56 
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This "natural" transition Reynolds number is too low for most situations, there- 
fore e is artificially set to zero until some other criterion is satisfied. A 

M 

Reynolds number on momentum thickness, Reg, is currently used to trigger transition. 
When a user prescribed Re^ is exceeded turbulent transport properties are in- 
troduced into the calculations; however they are reduced by a scale factor vary- 
ing between 0 and 1 to simulate a transition zone. The scaling factor, often 
referred to as an intermittency factor, has been reviewed in Reference 13 and a 
quadratic variation with streamwise coordinate has been recommended. However, 
due to the current state of transition data, a simple linear relation (Reference 14) 
was used . Thus 


e M = i(s) e M <ref) 


(40) 


where £,(ref) is calculated as if the riow were fuiiy turbulent and 
H 


1 (s) = — — 1.0 s. < s < 2 s. 

S t *■“ — ~ t 

1 (s) * 0 s i s t 

I (s) «= 1 s 2s t 


(41) 


2.3 COORDINATE TRANSFORMATIONS 

The equations of motion for a boundary layer flow can be solved in the 
physical (s,y) plane by numerous techniques, however it is generally advantageous 
to transform the problem to another coordinate system. The transformed coordinates 
offer the advantages of nondimensionalizing the solution, confining the solution 
to a narrower region, minimizing changes in the dependent variables, simplifying 
boundary conditions and occasionally result in the deletion of streamwise deriv- 
ative terms. This latter possibility occurs only under very restrictive sets 
of boundary conditions. The coordinate transformation in the present analysis 
is a variation of the Levy-Lees transformation and is derived in its entirety in 
Reference 1. The standard Levy-Lees transformation takes the form 
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C - / S Pl u l» 1 l r o ,C ds 


(42) 


r*u 


£ A* 


/2? 


The first alteration of this transformation is actually a mathematical conven- 
ience for carrying out the numerical solution. Introducing a stretching para- 
meter in the normal coordinate, a new coordinate system is defined by 


r- ? 


n-i 


(43) 


Ct H 


The parameter a is taken as a function of X only and is determined implicitly 
during the solution. Its purpose is to stretch the r) coo_3inate such that the 
boundary layer remains of constant thickness in the T) coordinates. 

Since a new variable a (£) is introduced, an additional relation is re- 

H 

quired. This is conveniently supplied by constraining seme arbitrary point near 
the boundary-layer edge, T) , to have a specified streamwise velocity, c, near 
(but something less than) the edge value: 


_ = cf ’ I 

In 


edge 


(44) 


where f is the transformed stream function defined as 


\ dn ■ “h r k 


dn 


(45) 
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and the prime denotes differentiation with respect to TJ so that 

< 46 > 

Examples of the utility of the stretching parameter a are contained in Reference 

H 

1. 

The second change in ths Levy-Lees transformation has to do with the 
transverse curvature effect. For very thin axisymmetric bodies, it is possible 
to have boundary layer thicknesses on the order of the body radius r Q . In this 
instance, it is necessary to treat r as a function of y, thereby including its 
variation through the boundary layer. The coordinate transformations become 

-5 - J Pl u l^l r o ds 

(47) 

r l f y ic . 

ft = — / pr dy 

« H ^ J ° 

Utilization of the above coordinate transformation relations results in 

A A 

a new set of governing equations in the (E,n) coordinate plane which will be given 
below. The hat ( A ) notation will be dropped for the remainder of the text for 
simplicity, however £ and r| are given by Equation (47) . Primes will refer to 
derivatives with respect to q expect when noted otherwide. 

The global continuity equation is automatically satisfied by the definition 
of a transformed stream function f(£,n), shown in Equation (45), and re-defined 
here in the final coordinate system: 


f 




(48) 
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0 V 
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p l u l w l r o 


d? 


The governing equations will be discussed separately. 
Streamwise momentum equation 



*= 2 



3f ' 

3 In C 


£' 2 


3 In a u 
& In £ 


3£ 

5 




) 


(50) 


In this equation, utilizing the technique of Reference 15, the transverse curva- 
ture effect is included entirely in the coordinate transformation and in the def- 
inition of t: 


t s 



1 + 


2a H /2T cosB 


u,r 


ro 




(51) 


where 0 is the angle between the surface normal and a plane normal to the body 
centerline (see Figure 2-1). Other definitions of interest are: 

C = -2iL (52) 

P 1 W 1 

9 In u. 

6 =' 2 imrr <53 > 

For solutions without consideration of transverse curvature, t is set to 1.0 
throughout the boundary layer. 
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Turbulent model equations 

The turbulent fluctuations are related to the mean field through the eddy 
models described in Reference 2. Eddy viscosity is described by a wall law and 
wake law, while eddy diffusivity and conductivity are related to eddy viscosity 
by turbulent Schmidt and Prandtl numbers: 


Defining 



(54) 

(55) 


6 5 


m. 


p l u l r c 


? - Pi 

5? 


RC 6 5 


P 1 U 1® 


i = 

M ~ 


(56) 

(57) 

(58) 

(59) 


The wall region eddy viscosity relation becomes 

(wall region) 


p(Re-) - 

e M - 2_ t *r 


p,a 


1 H 


? M * 0.018 




Re 


6 * 

X 


(wake region) 


(60) 


(61) 


where 


6 i " 6a H 



(62) 


Transverse curvature is not considered in determining the wake region length 
scale <5|. The governing equation for mixing length, which must be solved for 
the entire boundary layer although it is used only in the wall region, is 


l 
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(63) 
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Since mixing length is used only in the wall region, it is valid to use the one- 
dimensional expression for shear stress, Equation (24) . In transformed coordi- 
nates, this becomes 


T 

P 


^ fi Is. s_ + 

°H p L°H Re 6 


P W V W ft' 

p l u l 


(64) 


Energy equation 



3 *r . h» if 

3 In £ "T 3 In £ 


) 


(65) 


where q* is the normal iz.ed diffusive energy flux away from the surface including 
& 

turbulent fluxes 


q* - q a /a* 


( 66 ) 


The flux normalizing parameter a* is defined by 


a*. W£> 


(67) 


Diffusive energy flux q in the transformed coordinates is defined later in this 

£1 

section . 


"Elemental” species equations 



t 




3f \ 
J> In T ) 


( 68 ) 
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where 1* is the normalized diffusive flux of "element" k: 
J k 


J k ‘ 


169) 


Diffusive fluxes 


The normalized diffusive energy flux is given by 


q* » - Ji. 
* “h 


U 1 + pr T< 


+ C^RTy^ 


• + ±L Jc ♦ 

§c ^ l P v l v 2 ) 

3 + (h - h + c t RTy 3 )u^J 


r — • 

- ^ [*sf »i ♦ sf; T ' * si; <h ' - V' 1 ] 


(70) 


where Pr is the Prandtl number based on the frozen specific heat 


Pr 



(71) 


The turbulent contribution to the diffusive energy flux is contained in the last 
bracketed term, which is left uncombined with the other terms for clarity. The 
fact that the gross simplifications of the turbulent model are included in the 
same equation with the rather sophisticated unequal molecular diffusion model is 
merely a mathematical convenience stimulates by the requirement for calculations 
in all types of flow situations, including both laminar and turbulent flows. Un- 
equal molecular diffusion and thermal diffusion effects may be important in the 
laminar sublayer region of a turbulent boundary layer, however. 

Normalized molecular diffusive flux of species i is 


n 


_c 

a H §F 


\n 


+ a, 


- 


»] 


where Sc is a system property defined by 


SB = 


pDy 2 


(73; 


(74) 
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The Sc is a Schmidt numbe based on the self-dif fusion coefficient for a fictitious 
species representative of the system as a whole. The normalized molecular diffusive 
flux of the k th ' elemental" species is 


a «SH 




k T {Z k " K k )w 


»] 


(75) 


When certain groupings of parameters are constant so that the flow simi- 
larity assumption is valid, the terms on the right-hand side of the conservation 
equations (Equations 50, 65, and 68) vanish, in which case the conservation 
equations become ordinary differential equations. It should be emphasized that 
the equations as presented herein are equivalent to the corresponding boundary- 
layer equat'.ons presented in Section 2.1. That is, no similarity assumptions 
have been made in their development. 

Equations 67, 53, and 49 for a*, 8, and f , respectively, are indeterminant 
at the stagnation point of a blur'- body. Special forms for these equations valid 
at the stagnation point are shown in Reference 1 to be given by 


sp 


w. 


(^/ 6 V 

\ ' 'sp 


sp 


’ ( p wV a * ) 


(76) 


(77) 


sp 


where for Newtonian flow 


sp 


l/(iK + 1) 


(78) 


where i is the ratio of the crosswise to pitch-plane stagnation-point velocity 
gradients (i ** K for axisymmetric or planar flow) , and 


du l 

as~ 


(2P/p) 


'sp 


* /i 

sp / 


‘eff 


(79) 
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With R e £j an effective nose radius taking i ito account the shock shape, alter- 
natively, 6 and (du./ds) can be computed from curve fits of the inviscid 
sp J. sp 

pressure distribution. The transverse curvature parameter t also requires some 
special treatment at a stagnation point. The troublesome term is cos 0/r Q which 
is evaluated at a stagnation point by 



(80) 


In addition, to improve the accuracy of numerical integration procedures 
in the nose region, K and f can be computed by the following relations 


e 


1 

i +"1) 



d's 


2k+2j 


(81) 
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(2Q ~ >/2 
(k + 1) 


f V. <Ti) I 


rj 


which take advantage of the fact that u^/s and r Q /s x vary more nearly linearly in 
the stagnation region than do u^ and r Q . The basic approach is discussed more 
thoroughly in Reference 1 while the parameter i is discussed in Reference 16. 


2.4 BOUNDARY CONDITIONS 

The usual set of boundary conditions for the boundary layer flow problem 

consists of the specification of initial profiles for the dependent variables f', 

•>* 

H t , and K^, plus additional specifications of these quantities along the wall and 
at the edge of the boundary layer, and the specification of f along the wall. 
However, since the main utilization for the analytical technique presented here 
is to compute boundary layer properties for flows over ablating or transpired 
surfaces (heat shields, rocket nozzles, etc.), these boundary conditions have been 
greatly generalized. The numerous options resulting from this generalization are 
discussed below. 
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The boundary layer edge conditions typically are found from an isentropic 
expansion from known elemental gas composition and stagnation conditions. Thus, 
given a set of stagnation conditions and a description of loc?'. static pressure 
along the surface of interest, the techniques of Reference 4 may be used to es- 
tablish the entropy of the gaseous mixture which, when combined with the r- 
fied pressures, can be used to establish the complete equilibrium edge gas state 
at each body station. Edge boundary conditions then would consist of 


'■ 

f • 

edge 

““a 


* 

I. 

J 

H 

T 

* ^ 


j 

edge 

edge 

actual 

' 

** 

- K 


( 

edge 

edge 

actual 


where the subscript "edge" refers to conditions specified at , chosen to be 

outside the boundary layer (see Section 2.3). An additional constraint at the 
boundary layer edge which is necessary only when cubics are used is the require- 
ment o. zero slope, i.e.. 


edge 


0 


H' 

edge 


0 


(83) 


K = 0 

edge 

In addition to the specification of edge pressure, it is also possible to specify 
edge entropy and edge specie mass fractions to simulate the effects of entropy 
layer swallowing and nonequilibrium chemistry in the inviscid field. The tech- 
niques of Reference 4 are then used to establish the complete thermochemical gas 
state for nonisentropic , nonequilibrium expansions around a body of interest. 

Initial profiles of f", H^, and are more difficult to establish for the 
general problem, therefore calculations are often started with reasonable assumed 
profiles far upstream of the region of interest so that effects of erroneous 
assumptions will die out. Another possibility for initially laminar problems is 
to assume a similar solution as a starting profile. This assumption reduces the 
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equations to ordinary differential equations at the starting point, which may be 
solved simultaneously for a sat of profiles unique to the assumed edge and wall 
state. The similar solution is exact at a body stagnation point, therefore this 
option is particularly valuable for blunt body problems. 


The wall boundary conditions allow the widest selection of options. The 
simplest combination is the straightforward assignment of velocities, enthalpy, 
and elemental concentrations at the wall: 


f* = 0 
w 

f w - V» 

H T - h w «» 
w 


no slip 

specified P w v w 

specified enthalpy of gas 
at the wall 


(84) 






specified wall gas elemental 
composition* 


Hall to .peratures may be used to find wall enthalpy in the above formulation. 

Also, wall mass diffusive fluxes of up to three individual injectants may be 

assigned in lieu of and P w v w » With the values of the dependent variables 

w 

all directly assigned in this manner, the boundary layer problem is uncoupled 
from the surface chemistry interaction . 


The inclusion of surface material/boundary layer gas interaction chemistry 
in the boundary layer problem forms the second major set of wall boundary condi- 
tion options. Using the surface thermochemistry techniques of Reference 4, it 
is possible to specify given mass fluxes of the (up to) three injectants at the 
wall and require chemical equilibrium between the injectants, the wall material, 

and the adjacent gas stream. In this instance, the values of (i.e., T^) and 

~ V 

are found by simultaneous solution of the local surface chemical equilibrium 

w 

equations, surface mass balances, and the no-slip velocity boundary conditions. 
Alternatively, selected chemical reactions at the wall can be kinetically controlled 


It is physically unrealistic in most cases to^assign when diffusion coeffi- 
cients are unequal since the contribution to by w preferential diffusion of 
the various "elements" to the surface is not w known a priori. 



through Arrhenius-type rate law formulations (see sections 2.5 and 2.6) and in- 
cluded in the surface chemistry description. 

In the use of this boundary layer technique in conjunction with in-depth 
charring ablation analyses, the chemically active injectants might result from 
the pyrolysis of an internally decomposing material, surface material combustion 
or phase change, and mechanical icmovai. A variation of this type of wall bound- 
ary condition is to specify the wall temperature or enthalpy and allow the sur- 
face chemistry calculations to compute the necessary P w ^ w and 
the surface equilibrium wall boundary condition is 

f* = 0 no slip 

w 

specified P w v w 

(85) 

from surface equi- 
librium requirement 


The final wall boundary condition category involves the use of a steady 
state energy balance at the surface. A general surface energy balance can best 
be understood by examination of a schematic representation of the energy fluxes 
to an ablating or nonablating (m c «= 0) surface: 


f “ f (5) 
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Summing terms. 
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( 86 ) 


which is valid in either a transient or steady-state situation. In general, an 

in-depth charring ablation solution would be needed to provide the conduction 
• • 

term q cond and the pyrolysis gas rate, m^. Under steady state conditions, the 
internal pyrolysis "front" and the c-.itx.reu surface are assumed to be receding at 
the same rate, therefore requiring that the energy conducted into the wall mate- 
rial must equal the enthalpy rise of the wall material and pyrolysis gases. In 
equation form 


^cond 


- m (h - h ) - m (h 
c c c g g 
w ^ *w 


h C ) = 0 

g 


(87) 


Substituting into Equation (86) , the steady state energy balance becomes 
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In this equation, q is the wall value of the energy flux defined in Equation 

a 

w 

(70) , and is found in the course of the boundary layer solution. The surface 
equilibrium requirement is always used in conjunction with the steady state en- 
ergy balance. Therefore, if one specifies the compositions and heats of forma- 
tion of the pyrolysis gas and char materials, the simultaneous solution of the 
energy equation above and the surface chemistry relations mentioned earlier com- 
pletely couples the boundary layer flow to the surface response. The steady 
state assumption is good even in transient situations for large ablation rates 
or small thermal diffusivity of the ablation material (Reference 17) . In sum- 
mary, the use of the steady state energy balance results in the following; 


f' <= 0 no slip 

w 
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steady state 
energy balance 


( 89 ) 
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2.5 HOMOGENEOUS CHEMISTRY 

Chemical reaction rates are calculated using the procedures described in 
Reference 4. The m-th stoichiometric chemical reaction is written as 


I 


u. N . 
*3“ 3 



N . 
^3ra 3 


(90) 


and its reaction rate can be expressed generally by 




i a 


IT p> - kT TT p> 
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(91) 


The equilibrium constant can be determined from the standard state free energy 

change, AG°, to be 
Cl 


Ag 


In K = - 

p RT 
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(92) 


and the reaction rate constant can be expressed in the Arrhenius form 

A 


k = B T 
F m 
m 


\ r- -l 

m 

exp - E /RT 
m 


(93) 


The standard state free energy is a function of temperature only and is 
obtained for each molecular species from 


o o o 
Gj = Hj - T Sj 


(94) 


where enthalpies cure obtained relative to some chemical base state, often the 
elements in their most natural form at 298°K and one atmosphere (JANAF base 
state). If any other base state is consistently adopted, the Ag? will be unaffected. 
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The molar production rate of the specie i is then 


r. Ai? - y R \ R 
i Z-/ V im in ] m 
m ' • 

and the mass production rate is 


( 95 ) 




r i% 


(96) 


For the third body chemical reactions, the particular tnird body specie influences 

the production rate however in many cases, only the probability coefficient B is 

m 

different. In these cases, some reduction in the number of equations can be 

achieved by specifying a reference value of B and the relative efficiencies 

ra 

I* of each third body. Then by omitting the third body from the stoichiometric 
reaction an effective reaction rate k^, can be specified in terms of the molar 
concentrations n^ of the third bodies, i.e.. 



(97) 
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HETEROGENEOUS CHEMISTRY 


Heterogeneous chemical reactions are specified in a form analogous to 
homogeneous reactions noting that the surface, on which the reactions take place, 
is a third body. The calculation procedure is then similar to that for homogeneous 
chemistry. Of particular interest are surface catalyzed recombination reactions. 
There are two popular forms for specifying the effectiveness of a surface as a 
catalyst for Atomic recombination, namely, the catalytic efficiency and the 
catalycity k and is discussed in Section 4.2, The latter definition was chosen 

W 

because it could be related to the reaction rate k p and would therefore be con- 
sistent with the generalized surface chemistry calculation procedure. For in- 
stance, it is shown in the Appendix that for first order surface recombination 


k 

k = — 
F RT 


(98) 
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2.7 NODAL POINT DISTRIBUTION REFIT OPTION 

As will be shown in Section 5, nodal point distributions which are ap- 
plicable to equilibrium boundary layers are generally also applicable to non- 
equilibrium boundary layers. However, because o£ the increase number of “elements" 
that must be considered in nonequilibrium calculations and the attendently large 
matrix that must be inverted in the Newton-Raphson calculation procedure, sub- 
stantial reductions in machine time can be realized by reducing the number of 
•nodes- to a minimum. This is especially true when a lengthy laminar flow region 
is followed by transition to turbulent flow. In the laminar region, typically 
seven nodes are sufficient but in the turbulent region twelve to fourteen nodes are 
required. In addition a nodal distribution suitable for, say, the stagnation region 
may not be ideal for solutions far from the stagnation point . Thus a procedure 
was developed for adjusting the nodal point distribution and/or changing the number 
of nodal points to suit the local conditions where a solution is being sought. 


> 



The purpose of this procedure is to provide a means for maintaining an 
optimum nodal distribution for problems which include nonsimilar effects including 
transition to turbulence, blowing, entropi layer, pressure gradients, long stream- 
wise running lengths, etc. This readjustment is accomplished while preserving the 
fundamental characteristics of each profile, namely, basic profile shape, wall and 
edge derivatives, and integral properties. Potentially a number of bases may be 
identified for selecting nodal distributions and for making decisions relative 
to changing the existing distribution, e.g. mapping of any one of the velocity, 
temperature , and specie profiles. However, since adequate mapping of the velocity 
profile is the most commonly encountered problem, a selection criterion based on 
this parameter has been implemented, and the identification, evaluation and im- 
plementation of any other possible criteria has not been pursued at this time. 

Initially the selection criterion has been based upon maintaining a desired (specified) 
velocity ratio distribution across the layer; for nonsimilar turbulent flows then, 
for example, the nodal distribution, will change as a function of distance to account 
for the changes in velocity profile shape as the turbulent layor develops. The 
decision to refit is made following a converged solution and is based on whether 
or not the newly calculated velocities vary by more than a selected ratio from 
the desired values. 
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The REFIT procedure is currently valid for all allowed curve fitting 
options across the boundary layer, i.e. , all quadratics, quadratics with a final 
cubic and all cubics. It is also compatible with all of the entropy layer options. 
Finally as a result of the basic features of the REFIT option, it is possible 
to change the number of nodes used to describe the boundary layer. This latter 
capability has been programmed only for the case of transition from laminar to 
turbulent flow, as a means for eliminating the unnecessary and expensive extra 
nodes from laminar calculations. As such this option is limited to this application; 
however, potentially it may be programmed for more general application. The REFIT 
option is limited to a maximum of 15 nodes; however, as might be anticipated, the 
ability to maintain a more optimum distribution of nodes makes it possible to 
solve most problems using fewer nodes than normally required without REFIT. For 
example, for some long streamwise length, turbulent flows, it is either very dif- 
ficult or impossible to estimate in advance the best distribution for the entire 
length using all 15 nodes. With REFIT, it is possible to achieve good results 
with minimal selection of desired velocity ratios using 12 nodes. Since solution 
times vary roughly as the number of nodes squared, this represents a saving of 
40 percent in computer time, some of which is used in the refitting operation. 
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SECTION 3 

INTEGRAL MATRIX SOLUTION PROCEDURE 

The solution of the transformed boundary layer equations presented in 
Section 2 uses an integral matrix method which has been developed specifically 
for the solution of chemically reacting, nonsimilar, coupled boundary layers. A 
complete presentation of the integral matrix procedure was included in Reference 
1, where solution of laminar flow problems was discussed. In the present effort, 
this technique has remained essentially unchanged, however new variables and 
equations have been added to describe the nonequilibrium aspects of the flow. 

The present discussion will therefore review only the highlights of the method, 
and the reader may refer to Reference 1 for more details. 

In the integral matrix procedure, the primary dependent variables and their 
derivatives with respect to n are related by Taylor series expansions such that 
these dependent variables are represented by connected quadratics or cubics 
(either option is available). That is, f', H^, and are expanded in Taylor 
series form and the series are truncated to reflect the proper polynomial repre- 
sentation. A nodal network is defined through the boundary layer and the Taylor 
series expansions are assumed valid between each set of nodes, with an additional 
requirement of continuous first and second derivatives (a spline fit) . Primarily 
for convenience, the conservation equations are integrated across each "strip" 
(between nodal points) using a unity weighting function. The linear Taylor 
series expansions together with linear boundary conditions form a very sparse 
matrix which has to be inverted only once for a given problem. The nonlinear 
boundary layer equations and nonlinear boundary conditions are then linearized, 
the errors being driven to zero using Newton-Raphson iteration. 

3.1 integral strip equations with splined interpolation functions 

Consider the boundary layer in the region of a given streamwj.se station s 
as being divided into N-l strips connecting N nodal points. These nodal points 
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are designated by ru where i = 1 at the wall and H at the edge of the boundary 
layer. Consider a function p(n) which with all its derivatives is continuous 
in the neighborhood of the point rj = n^. Then, for any value of T) in this neigh- 
borhood, p(n) may be expressed in a Taylor series expansion as 


p i + l 


P* + P !6n + p'| 


+ p r l«nLl 

2 6 


+ p 1 !" 


(6n) - . 
~2l 


where 


«n n 1+1 - 


( 100 ) 


Conventional finite difference schemes, in effect, typically truncate the 
Taylor series after the first term and use the resulting expression to relate p' 
to p, etc., that is 


P i+1 ~ p i 

6n 


( 101 ) 


Round-off error is then of order (fin) 2 and many nodes must be chosen to bring 
this value down to acceptable limits. One can achieve a reduction in the number 
of nodes for a given accuracy by employing a quadratic or cubic relation repre- 
senting the function p over the interval of interest. This can be achieved by 
truncating the Taylor series after the third or fourth term. The cubic approxi- 
mation will be used for the remainder of this discussion. The p^ can be consid- 
ered to be any of f., f^, fV, f£’, H T , H.J, , , or K£ . Since the highest 

derivatives of the dependent variable! which appear in the boundary layer 
equations are fV, Hjj, and K£ , it is reasonable to truncate the series at the 
next highest derivatH 

T). and T). , , , that is, 

1 1+1 


next highest derivative and to consider that derivative as being constant between 


i f iVi 


f iVi - f r 
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Thus, rather than using finite difference approximations similar to Equation 
(101) which are substituted directly into the governing differential equations, 
a set of 1 inear relations between the dependent variables and their derivatives 
is obtained and is solved simultaneously with the governing differential equations. 
These linear relations are of the form 


- f i+i + f i + n 6 " + + f i 


». 8 ... (6n) * 

8 + f i+i TT~ 


= o 


' Pi+i + Pi + Pi 6r > + Pi 


•6n + p? IMl + p ^ +1 Ji-5l 


- Pi.l + Pi + Pi Q * Pi 


i+1 T 


in = o 


(103) 

(104) 

(105) 
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\ 
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where in Equations (104) and (105) the p^ represents f|, 
sets of . 



and each of the K 


i 

Notice that f' has been taken to be a cubic over each strip, rather them 
the stream function, f, since it was desi-ed to represent velocity (u = u^'/Ojj) 
with the cubic. Equations (103) through (105) above, when written for each 
adjacent pair of nodes, give (3 + 2K) (N - 1) simultemeous algebraic equations 
for the N (4 + 3K) + 1 unknowns, f , f', f", f’", a , H , H’ , H" , K. , K.' , 

K£ at each streamwise station, where K is the number o? elementa? species.* 

The Taylor series equations are written for only K-l species since the overall 
mass balance equation supplies the remaining elemental concentration. Additional 
relations- must come from the governing differential equations and the boundary 
conditions. It is import 3m t to note that the f, f’, etc., are treated as indi- 
vidual variables related by algebraic equations. It is also important to note 
that the coefficients in Equations (103) through (105) are functions of 6r) only? 
therefore, this portion of the resulting matrix need be inverted only once for a 
given problem. 


The mixing length is not included in this variables count since mixing length 

(as well as e„ in the wake region) is treated as a state property. 

M 
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The conservation Equations (50), (65), and (68) contain streamwise deriva- 
tive or "nonsimilar" terms. In the present solution technique, two or three 
point fin' te difference formulas are considered sufficient to express these de- 
rivatives, since gradients in this direction are not severe. As in Reference 1 

2 [ Uln 4 * V >t * d l ( >t-l * d 2< >t-2 < 106 > 


where ( refers tc the previous streamwise station, 


d 


£“£- 1 


* d. 


A- 1 


' d 2 " 0 


for two-point difference and 


, , A-i + A - 2 

a. = l " V 
0 £*£-1 £ u £-2 




A- 2 


£°£-l £-1 q £-2 ' 


d„=2 


A-l 


2 A-2 £-l A £-2 


(107) 


(108) 


for three-point difference where typically 

A-l = ln " ln ^JL-1 = ^AA-l* <109) 

The three-point difference relation is generally used unless a similar solution 
is desired (in which case d o = d^ = d£ = 0) or unless the point in question is 
the first point after either (1) a similar solution or (2) a discontinuity 
(e.g., where the body changes shape abruptly, or wl:are mass injection is suddenly 
terminated) . 

The next step in the treatment of the conservation equations is their in- 
tegration across the boundary layer "strips". The primary reason for this inte- 
gration is to simplify the r)-derivative terms in the energy and species conser- 
vation equations, since it is not convenient to express the complex q* and j* 
terms in derivative form. The solution can actually proceed very nicely with- 
out integrating across, strips (see Reference 8) without any noticeable change 
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in speed, accuracy, or stability for simplified problems such as incompressible, 
nonreacting flows. The weighting function for integration between nodes in this 
integral method is unity. In the terminology of the general method of inte< ~al 
relations, where integrals are carried form 0 to 00 in r| (Reference 18) , a rquare 
wave weighting function is used which is unity across the strip in question and 
zero elsewhere. The equations are then integrated N-l times with the square 
wave applied to each strip in succession. Using the momentum equation as an 
example, the integration from i-1 to i results in 



f ’ 2 dn 


i i 

- f £•(«„£’ + Vi-i * Vi-*''”' - f *-[V» “a 

•'i-l *'i-l 
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+ d x (ln « H ) + ^2 dn cl H ) £-2] dn “ f f " (d o f + a l f /.-] 
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+ d 2 f £ _ 2 )dn 


( 110 ; 


The T&ylor series approximations introduced earlier can also be used to express 
the integral terms above. ?.s demonstrated in Reference 1, the term J* 1 f’p di) 
becomes * 




f‘ p dq 
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where 
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This technique is used to rewrite each of the integral terms in Equation (11G) 
above of the formj* 1 f'p dq. The remaining integral term in th*. momentum 
equation, J* 1 (p^/jlijclq is evaluated by approximating these functions as cubics 
over the strap and integrating directly. This yields 



p l p i-l\ 6n 2 
.2 I 12 


(113) 


The production term is assumed to vary linearly across the strip so that the in- 
tegral of <t K /p is 



(114) 


These approximations are not quite as good as the approximations for f', I1 T and 
since continuity of derivatives is not guaranteed at the nodal point. 

Direct substitution of these approximations for integral terms into the 
governing equations results in the following forms. 
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Energy 
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.The following definitions are necessary: 
ZP 
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with 


**1 C d l p l-l.i + d 2 p l-2,i 
W 2 “ d l p t-l,i + d 2 p t-2,i 
W 3 “ d l p I- 1#i + d 2 p l-2,i 
W 4 = d l p l-l,i-l + d 2 p l-2,i-l 


and is defined adjacent to the brackets in each term that uses these definitions. 

•Hie conservation equations provide (K+l) (N-l) more equations for the 
N(3K + 4) + 1 unknowns, thereby closing the problem. However, before discussing 
how this set of algebraic equations is solved.- Section 3.2 describes in detail how 
the mixing length differential equation is solved. 

3.2 SOLUTION OF THE MIXING LENGTH EQUATION 

The mixing length equation is a first order linear differential equation 
whose solution can be written directly in general terms. The differential 
equation is 


dl 
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(Xa H n - 1 ) 


(63) 


Defining 
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( 120 ) 


results in 
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The solution to this equation is 


1 — Kci 
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The remaining problem is to evaluate the integral terms. Defining 


L(n) = ^ 


. _ f n P d"’" 

/ '• J o 

A e d ** 


r n p dn* 


(123) 


yields 


A = KOtjj (n - L) 


024) 


Reference 2 presents a complete description of the technique used to evaluate 
L(r?). In essence.. P(n) is assumed to vary linearly over the interval n. , to 

X“i 

and the integrals are expressed in a more tractable form. The final expression 


is 



(125) 

(126) 

(127) 

(128) 

(129) 
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The Dawson Integral, D w ( ) , can be evaluated from tables (Reference 13) or by a 


• • if •> ^voxua^xun 


'thod is used ir. ths present snsiysis . *ry>vi<s 
combining Equations (124) and (125) , an explicit recursion formula for mixing 
length at each node is obtained. This mixing length is a function of local shear, 
viscosity, and density through the variation of P( ), and is re-evaluated at each 
node on each iteration during the course of a solution. 


3.3 NEWTON-RAPHSON ITERATION FOR A SOLUTION 

A complete description of the Newton-Raphson iteration procedure as ap- 
plied to the laminar equations of motion was given in Reference 1. Since the 
procedure is basically unchanged with the addition of nonequilibrium chemistry 
it will be reviewed only briefly here, with emphasis on the recent additions. 

To illustrate the Newton-Raphson method, consider two simultaneous non- 
linear algebraic equations 


?(»,«) a= 0 fi(x.v) = 0 


(130) 


the solution for which is given by x = x, y = y. Define x and y as the values 

th “ _ m 

of x and y for the m iteration. The desired solution f(x,y) can be expressed 

in a Taylor series expansion 
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3G < x m' y m 5 
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(131) 


The Newton-Raphson method consists of replacing (x,y) by (x y , ) on the ■ 

mrrX nH-x 

hand side of these expressions and neglecting nonlinear terms in and 
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y m+ ^ - v^. This yields the set of simultaneous equations 
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( 132 ) 


U33) 


Ax = X - - X 
m m+1 n 


K = F „*l - y m 


(134) 


The Ax and Ay are the corrections to be added to x and y , respectively, to 
mm m th m 

yield the values of the dependent variables for the m + x iteration. Hera 

F(x ,y ) and G(x ,y ) are values of the original functions F(x,y) and G(x,y) 
m m m m 

evaluated for x = x and y = y . As the corrections approach zero, the F(x ,y ) 

mm mm 

and G(x m ,y m ) approach zero. Hence, it is appropriate to look upon these as 
errors associated with the original Equation (130) . It is apparent that this 
procedure can be extended to an arbitrary number of functions and a corresponding 
number of primary variables. 

For the purpose of the present analysis, it has been Found most convenient 
to consider the primary variables as f^, f£, fV, f?’, , H.J, , , 

K£ , and This amounts to (3K + 4)N + 1 unknowns where N is the numfeer o£ 

noctes and K is the number of elemental species to be considered in the boundary 
layer. Recounting the number of equations, we have 
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Eqn. Numbers 


No. of Equations 


[•»y mj ssriss expansions (Iim) - (iub) (N — 1) [5 + 2(k ■ iil 


Boundary layer equations (115) - (117) 

Boundary conditions (82) , (83) , (84 


(82), (83), (84) 
or equivalent 


definition 


Total 


(N - 1) (K + 1) 
3K + 4 


N(3K + 4) + 1 


Other secondary variables such as e, p, T, etc. are expressed in terns of those 
listed above. The corrections in these secondary variables are therefore found 
in terns of the corrections to the primary variables. 

The use of the Newton-Raphson technique for the current set of equations 
requires the evaluation of the partial derivatives of each equation with respect 
to each variable. The partial derivatives of the Taylor series equations and 
linear boundary conditions are exactly the same as in Reference 1. The deriva- 
tives of the conservation equations are: 

Momentum 
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where the ERROR is given by the left-hand side of Equation (13 5) evaluated for 
th . 

m iteration. 
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where the ERROR is given by the left-hand side of Equation (11 fij for the hi 
iteration and Aq* is given by 
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where the ERROR is given by the left-hand side of Equation (117) evaluated for 
the m iteration and Aj* *" s given by 
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The technique of relating corrections on secondary variables such as C, 
o, T, Pr, etc., to corrections in primary variables was fully explained in 
Reference 1. The same techniques are used for the corrections At and Ae anu 
M^p). 

Once the correction coefficients (partial derivatives v,ith respect to 
each primary variable) for each equation at each nodal point are found.- they 
are arranged in matrix form for further manipulation. The order of the primary 
variables and the order of the equations is of some importance in the matrix for- 
mulation. It is most convenient to divide the variables into "linear" (symbol L> 
and "nonlinear" (symbol NL) sets, namely 
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where the linear equations are the Taylor series equations and some of the bound- 
ary conditions. The purpose of the partitioning is to allow operations on sec- 
tions of the coefficient matrix which result in significant simplification of 
the overall inversion. In particular, since the coefficients of the linear equa- 
tions are all constant or functions of the fixed nodal spacing, this portion of 
the matrix (the AL portion) can be diagonalized once and for all in any given 
problem. In essence, the corrections on the linear variables AVL are always ex- 
pressed in terms of the nonlinear variable corrections Avnl. The choice of 
linear and nonlinear labels for the variables is somewhat arbitrary, but care 
must be taken that the AL matrix not be singular. It has been found convenient 
to arrange the variables into the linear and nonlinear groups as follows: 

L H (iH T ' 

Shen arranged 
); Avnl h (Ah^, , 

Ah t » Ah t , . . . Ah t ) ? and K-l sets of Avnl r (Ak£ , Ak^ , Ai^ , . . . Ai^ ) . w 
The W order 2 of the linear equations in the present matrix procedure is: n ^ 
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Ah£ ); and K-l sets of AVL^Al^ 


Ak^ , . . . Ak^ , Ak» , Ak^ , . . . Ak£ ) . The nonlinear variables are 
in the following order: 2 Avnl n (Aa , Af , Af", Af ' , Af!,... Af' 


No. of 

Equations Description of Equations 

3N-2 Linear Boundary conditions and 

Taylor series for f, f', f", f'” 


2N Linear boundary conditions and 

Taylor series for H^, H^,H^ 

(K - 1) (2N) Linear boundary conditions and 

Taylor series for K^, K£, K£ 


The nonlinear equations are sequenced as follows: 


No. of 
Equations 

4 


N - 1 


Description of Equations 

Nonlinear boundary conditions 

and a constraint 
H 

Momentum equation for each pair 
of ncdes 


N Energy equation for each pair of 

nodes plus wall enthalpy equation 

(K - 1) (N) K-l sets of "elemental" species 

equations for each pair of nodes 
plus wall species equation 
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Special logic has been written for the matrix inversion, taking advantage 
of the regular sparseness of the matrix. Once the corrections for the linear 
and nonlinear variables are found, these corrections are added to the variables 
to form the new guesses. The magnitude of the errors for each equation are 
checked and the procedure advances to the next iteration if the absolute values 
of the errors exceed prescribed upper limits. If the errors are acceptable, 
iteration is completed for the current streamwise position £. Typically, three 
to six iterations are required to reach a satisfactory solution. 


SECTION 4 


CONSIDERATIONS FOR SHUTTLE APPLICATIONS 
4.1 TRANSITION TO TURBULENT FLOW 

For the shuttle vehicle, transition influences not only the choice of TPS 
material but also sizing of the TPS. Because transition criteria are not well 
established, many studies have been conducted to assess the affect of various 
transition onset criteria on peak temperatures and TPS weights. For instance, 
Reference 20 considered the RI 134B vehicle for both constant and variable angle- 
of attack entry trajectories and found that 95 to 99 percent of the TPS weight 
is determined by laminar flow heating. The remainder is due to transition and 
turbulent heating. This 1 to 5 percent requirement for turbulent, flew appears 
small, however for an P f ’I system, the total TPS weight may be 20,000 pounds re- 
sulting in a 1,000 pound variance depending on when turbulent flow will occur. 

The reentry heating conditions determine the peak temperatures experienced 
by an equilibrium radiation surface which in turn influences the width of the 
entry flight corridor (References 21-24) (along with equilibrium glide and maximum 
force constraints) . Not only are the conditions which influence transition onset 
not well established, but the length of the transition zone is not clearly defined. 
Current best estimates place the ratio Re tj ,/Re t at between 2 and 3. Thus a large 
portion of the vehicle may be in neither laminar nor turbulent heating but 
rather an in-limbo transitional state. 

The question of what constitutes transition has been discussed by several 
authors (for example. References i3, 25-28). Clearly a general knowledge of turbu- 
lence and the transition from laminar to turbulence is not available. However 
the intermittancy concept, which envisions local flow conditions which are inter- 
na ttantly laminar and turbulent (References 13, 26, 29, 30) (an average condition 
being the transition state) has received favorable analytical attention. In effect, an 
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intermittancy factor is employed which scales the turbulent eddy viscosity between 
0.0 and 1.0 of its fully developed value. Both linear as well as nonlinear 
dependences on flow length from the onset of transition have been considered. 

Many possible variables or parameters can affect transition onset. The 
fact that surface roughness, gaps, discontinuities or external trips can influence 
transition onset is classical; also heat transfer to surfaces have been known to 
have a stabilizing effect on the boundary layer. Other variables which have been 
shown to induce transition include nose bluntness (Ref. 25, 29, 31, 33) masi injection 
(Ref. 29, 33), entropy swallowing (.Ref. 25), adverse pressure gradients (Ref. 2'', 34) 
boundary layer edge conditions (Ref. 25, 34) free stream unit Reynolv s numbers (Ref. 25, 
27, 33), and vehicle angle-of-attack (Ref. 25, 33). In addition, iu wind tunnel ex- 
periments, tunnel size, wall effects and noise (Ref. 33) must also be considered. 

Because of the large number of variables that apparently hrva influence 
on transition onset, a fundamental theory which encompasses all of these variables 
is not available. Many correlations of flight and wind tunnel data have been 
attempted and for shuttle applications, the correlations of References 25 and 35 
appear to agree that the boundary layer edge conditions, free stream unit Reynolds 
number and momentum thickness are good correlation parameters. Similar conclusions 
were obtained in the Philco-Ford (Ref. 36) correlations presented in Reference 29. 
Reference 29 showed that both the Philco-Ford and the McDonnell-Douglas (Ref. 35) 
correlations follow the trend of a large number of flight data. However, most of 
th c data was for slender cones with ablating nosetips and a separate, more involved 
correlation was obtained by Martellucci (Ref. 29) which includes the effects of nose 
bluntness and mass addition. For shuttle applications, the amount of mass addition 
will be small; in addition no significant shape change is expected so that the simpler 
correlations of References 25, 35 and 36 are recommended. 

Reference 25 presents a correlation of thfc form 

Re. ■ f (M /Re/ft‘) 141 

t e 

Where Re is the transition Reynolds number based on edge properties and running 

length from the stagnation point, M is the edge Mach number and Re/ft is the local 

6 

unit Reynolds number also based on edge properties. This correlation was obtained 
for a simulated shuttle configuration under wind tunnel conditions (M^ « 10, 

1.0 X 10® <_ Re^/ft <2.4 X 10 6 ) for angles-of-attack between 5 and 35°. From the 
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data in Reference 25 the functional relationship (141) can be approximated by the 
expression 


Re t » (M^/Re^/ft ) 1,85 (142) 

Note, however that the above correlation was obtained for a prescribed vehicle 
configuration at ® 1 ' for a particular wind tunnel. 

The correlations of Reference 35 are for a number of delta wing configura- 
tions under wind tunnel test conditions and has also agv ed with slender cone flight 
data. This correlation has the form 

Re- « f(a)(Re/-)°* 2 h (143) 

e t 

Where f(a) is a function of the vehicle angle of attack, Reg is the Reynolds 
number based on momentum thickness, is the edge Mach number and Re/x is the 
local unit Reynolds number based on dge properties. For the data used to obtain 
the correlation (143), f(a) has a value of ~ 10.0 for a < 35® and increases to 
" 6 for a ■ 60®. The accuracy of the correlation, as indicated by data scatter, 
is about a factor of 4. If the unit Reynolds nomber effect xs not considered, 
the same data is represented by the correlation 

Re 0 fc * 150 M fe ( 144 ) 

with data scatter as high as a factor of 6 . 

The dependa :e on unit Reynolds number is small (to the 0.2 power) and 
Martellucci's independent comparison witn t’e same and other data show a larger 
amount of data scatter which led Martellucci to question vhether o’* not there 
is a unit Reynolds number dependence. In view of the fact that Reference 25 and 
35 both observed a unit Reynolds number effect for shuttle type vehicles, either 
correlations (142) or (144) should be used as transition criteria bearing in 
mind however that the data scatter has an uncertainty -actor of about 4 and 
that the unit Reynolds number effect, if real, is apparently small. Finally it 
is suggested that some form of intermittancy factor be incorporated to simulate 
a transition zone and that the length of this zone be determined by Re T /Re t ~ 2 . 0 . 


4-3 




l 


l 


i 


i 


4.2 SURFACE CATALYZED REACTIONS 

Since significant dissociation will occur in portions of the flow field 
surrounding a shuttle vehicle, surface catalyzed recombination reactions repre- 
sent a significant portion of the energy transfer to the vehicle. The reduced 
heating which can be realized from low catalycity surfaces would locally reduce 
the equilibrium radiation temperature for any given flight condition. However, 
if a rapid change in catalycity from very low to very high occurs as the flow 
progresses downstream, then the downstream heating is increased by virtue of the 
noncatalytic upstream section. It is clear that optimization would require a 
trade-off study which can be performed with the current code provided that adequate 
data is available for catalytic efficiencies of materials of interest. Surface 
catalyzed reactions are often defined in terms of catalytic efficiency y or cata- 
lycity K w . 

The catalytic efficiency is defined as 


Y = 



<145) 


Where N is the collision rate of atoms with the surface and N p is the rate at 
which these atoms reccmbine due to collision with the surface. If heterogeneous 
surface catalysis were due to the simultaneous collision of two atoms and the sur- 
face, the recombination rates would be similar in magnitude to homogeneous reac- 
tion rates. However, very large catalytic rates (y ~ 1) have keen observed which 
lead to the postulation of mechanisms which require an adsorbed layer of atoms 
on surface active sites with reaction formulas that are identical to homogeneous 
reaction formulas (Ref. 37-39). Designating A* as a surface activation site, 
then the Langmuir-Hinshelwood mechanism is 


C + A* -P* OA* 

OA* + OA* -Pf 0 2 + 2A* 


Langmuir-Hinshelwood 

(L-H) 


and the Lideal-Eley mechanism is 


OtA*^ OA* 

0 + OA* 0 2 + A* 


Rideal-Eley 

(R-E) 
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In all likelyhood both mechanisms occur simultaneously with rates that depend 
on the availability of surface active sites and the mobility of absorbed atoms. 

At high surface temperatures, surface diffusion rates are high and the L-H mech- 
anism is expected to predominate whereas at low surface temperatures the direct 
collision R-E mechanism is expected to predominate. In any event, it is possible 
for both mechanisms to be rate limited by the reaction 

0 + A* OA* | 

« 

which depends on the concentration and state of atoms near the surface since the 
collision rate is a function of the partial pressure of the atoms and their tem- 
perature and the probability of the atoms adhering to the surface depends on the 
kinetic energy (or temperature) of the atoms. Both mechanisms, in terms of the 
concentration of atomic species, cam exhibit first or second order behaviors. 

For instance, in the L-H mechanism, a slow adsorption of 0 atoms and a rapid 
surface migration of OA* will result in a first order reaction; similarly, the 
R-E mechanism will be first order for a rapid adsorption of 0 atoms which forms a 
high surface density of OA* such that the reaction 0 + OA* 0^ + A* controls. 

The catalycity k^ is defined (Ref, 40) for convenience in gas dynamic studies 
in terms of local atom concentration and diffusive flux. For a cata?v<-ic reaction 
the diffusion rate, j , is equal to the surface reaction rate so that the net 
overall reaction for both the L-H and R-E mechanisms can be written as 


or 


20 (Second order) 

0 1/2 0 2 (First order) 


(146a) 
C 46b) 


For nonequilibrium boundary layers, the forward direction is the most probable. 
If the heterogeneous reaction were treated as a homogeneous reaction (with the 
surface as a third body) then the reaction rate (146a) has the form 


(M 0 > = - k f [(M 0 ) 2 - )] 


(146c) 


Where (M^) is the molar concentration of specie i, (MJ is its production rate 

per unit surface area, and k is the molar equilibrium constant. For low surface 

c 

temperatures* and a nonequilibrium condition near the surface the term (M 02 ) A c will 

be small compared to (M ) 2 so that the reaction rate can be approximated as 

o 


Low temperatures in the current context is defined as a condition in which dissoc.- 
ation is negligible. 


j 


! 

t 

4 
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(146d) 


(M ) = -k (M ) 2 
o f o 


For near-equilibrium conditions the reaction rate will approach zero with (M ) 5 

o 

- (M )A and can be achieved with appropriately large values of k . Thus 
02 c r 

at low temperatures reaction (146d) would still be valid since (M ) approaches 

o 

zero. 


The catalycity k^ is defined for convenience in terms of the local atom 
concentration and diffusive flux. That is, for a catalytic reaction the 
diffusion rate, j^, is equal to the surface reaction rate. Then making use of 
reaction (146d) , with an appropriate transformation to mass production rates, we 
have (Reference 40) 


‘Vw’X ' J. 



(147) 


Where n is the order of the reaction and c is the mass fraction of dissociated 
specie i. 


Using a kinetic theory definition of the collision rate, i.e. , 


N 


(2mn o feT) 




Where the subscript o is for dissociated species, y and k^ are then related by 


/ 2nfc y /i / * \ n> ~ 1 . 

( o \ [ o \ f »n-l 

w v KT W / ) p ° 


(148) 


The definition of k is convenient since, with it, the surface reaction rate can 

w 

be expressed in th? acceptable form for a forward reaction, i.e., reactions (146a) 
or (146b). 


Note that, for first order reactions (n ■» 1) the relationship between y 
and k w is independent of the specie partial pressures. However for n > 1, the 
partial pressure is required in linking y to k^. in addition, if the reaction 
is first order, at any given value of T^, there is a maximum finite value for k^ 
which corresponds to y * 1. There is then an obvious discrepancy in t l e. relation- 
ship (148) since, as a reaction rate, the limits of k^ should be between 0 and 00 . 
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4. 3 HOMOGENEOUS KINETICS 

Very comprehensive compilations of kinetic reaction rates for oxygen- 
ni trog en-carbon-hydrogen systems are presented in References 41 through 44. Some 
of these rates as well as rates from other sources are shown in Table 4-1 and 
4-2. For shuttle applications the oxygen-nitrogen reactions are most important; 
the introduction of small quantities of carbon and silicon compounds , due to 
ablation or surface oxidations would have only a small effect on the boundary 
layer solution. Thus oxygen-nitrogen reactions which are important for shuttle 
environments are shown in the first table and selected carbon-oxygen-nitrogen 
reactions are shown in the second table. All rates are presented in the modi- 
fied Arrhenius form 

k f - AT* exp£-|i| (149) 

with T in °K and the units of k^ consistent with (moles - cc - sec) . Only for- 
ward rates are shown; the presumption being that reverse rates can be calculated 

from the equilibrium constant k or k which in turn can he calculated from 

pc 

free energy considerations. 

As noted by Dryer (Ref. 45) Equation (149) is not always the best correlation 
of data, especially over a wide temperature range. T> is possibly one of the 
reasons why Reference 43 recomnends two reaction rates; ne applicable at high 
temperatures where the endothermic reaction is dominant and one applicable when 
the exothermic reaction is dominant. Although it is possible for both exo- and 
endothermic reactions to be important in different parts of the flew field, only 
the high temperature values of Reference 43 are shown in Table 4-1. 

At wall temperatures of interest, carbon sublimation will be negligible 
and the primary ablation product from a carbon heat shield will be CO. The 
relative reactivity of CO with air species was compared and Table 4-2 includes 
only those reactions which are most likely to havp a significant influence on the 
boundary layer solution. 
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TABLE 4-1 
AIR REACTIONS 







A 



N 

E 

REF. 


l) 

°2 

+ 

M ^ 0 + 0 + M 





(°K) 






M»°2 

3.3 

(19) 

- 

1.0 

59,400 

43 





M = 0 

9.0 

(19) 

- 

1.0 

59,400 

43 





m = n 2 

7.2 

(18) 

- 

1.0 

59,400 

43 





M = N, NO, Ar 

3.6 

(18) 

- 

1.0 

59,400 

43 





M = N, NO, Ar 

2.5 

(16) 

- 

0.5 

59,400 

46 

> 




M * N, NO, Ar 

3.5 

(18) 

- 

l.C 

59,400 

47 

1 

2) 

N 2 

+ 

M ^ N + N + M 











M = N 

4.15 

(22) 

- 

1.5 

113,100 

43 





M « N 2 

4.75 

(17) 

- 

0.5 

113,100 

43 





M = 0 2 ,0, NO, Ar 

1.92 

(17) 

- 

0.5 

113,100 

43 





M = 0 2 ,0, NO, Ar 

2.0 

(20) 

- 

1.5 

112,500 

48 





M = GENERAL 

4.75 

(17) 

- 

0.5 

112,500 

46 

£ 

3) 

NO 

+ 

M ^ K + 0 + M 







* 




M - NO, 0, Ar 

7.8 

(21) 

- 

1.5 

75,600 

43 

t. 




M = 0„, N 
2 2 

4.0 

(20) 

- 

1.5 

75,600 

43 





M = O a , N 2 

5.5 

(20) 

- 

1.5 

75,000 

48 





M - 0 2 , N 2 

5.3 

(16) 

- 

2.5 

75,000 

46 

; 

4) 

NO 

+ 

0 +? N + 0 2 







■< 





3.2 

(9) 


1.0 

19,550 

46, 48-51 

- 





3.2 

(9) 


1.0 

19,700 

43 

\ 

5) 

N 2 

+• 

0 ^ NO + N 

6.6 

(13) 


0.0 

37,750 

52 -55 

* 

F 





7.0 

(13) 


0.0 

37, ”50 

46, 50, 51 

S 





6.75(131 


0.0 

37,500 

43 






6.8 

(13) 


.0.0 

37,500 

48 

i 

l 

6) 

N 2 

+ 

0 2 NO 







£ 





9.1 

(23) 

- 

2.5 

64,250 

48 






8.4 

(13) 

- 

0.5 

61,600 

43 


i 
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TABLE 4-2 


1) C0 +M<-*C0+0 + M 


2) CO + 0 +* C + 0 2 


3) CO + 0 2 *.-* C0 2 + 0 


4) CO + N C + NO 


5) CO + N 2 -e* CK + NO 

6) CO + N CN 4- O 


7) C + NO CN + O 


8) C + N 2 «-■* CN + N 


9) C + O-tM^CO + M 


AIR AND CARBON REACTIONS 


A 

N 

E 

REF. 



(°K) 


2.3 (13) 

0.0 

37 ,050 

56 

1.8 (20) 

- 0.7 

65,000 

43 

1.3 (18) 

- 0.58 

62,290 

48 

3.0 (12) 

0.5 

132 

43 

2.5 (12) 

0.0 

48 

57 

1.6 (13) 

0.0 

41 

58 

7.0 (12) 

0.0 

51 

49 

3.6 (12) 

0.0 

50 

59 

4.8 (17) 

1.0 

100 

43 

6.0 (11) 

0.0 

276 

43 

2.4 (12) 

0.5 

92 

43 

6.0 (11) 

0.0 

0.0 

43 

1.2 (12) 

0.0 

56.8 

43 

3.0 (16) 

- 0.5 

0.0 

60 


f 
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SECTION 5 

TYPICAL RESULTS BASED ON EQUILIBRIUM EDGE BOUNDARY CONDITIONS 


Because of the difficulties of attaining inviscid f lev: solutions with 
nonequilibrium chemistry for general bodies, many boundary layer calculations 
are performed using equilibrium edge conditions and, in the case of blunt 
bodies, isentropic expansions from the stagnation point. For shuttle vehicles, 
which decelerate at high altitudes and have long characteristic lengths compared 
to their nose radix, these assumptions are not valid. Nevertheless, thes-* 
assumptions are convenient and were used to establish the significance of 
variables such as nodal point distribution, surface catalycity, homogeneous 
chemistry, and pressure. In most cases to be discussed stagnation point solutions 
on a nose radius of 2.325', a total enthalpy of 9013 Btu/lbm and stagnation 
pressures between .004 and .09 atmospheres were used since this would be repre- 
sentative of a shuttle vehicle at relatively high altitudes. The parameters of 
all cases in this section are shown in Table 5-1 and Figure 5-1. The homogeneous 
reaction rates are shown in Table 5-2. The values shown in Column 2 of this table 
were used in most cases in this section, however, column 1 was used for com- 
parisons shown in Section 5-1. 


5.1 EFFFCT OF HOMOGENEOUS KINETICS 

The effect of two different sets of kinetic data (Table 5-2) are shown 
in Figure 5-2 for noncatalytic walls. A significant difference is noted in the 
N concentratic and results in an equally significant difference in heat 
transfer ( -v 15% from Table 5-1) . Although not shown, the difference decreases 
as the catalycity increases with no significant differences for fully catalytic 
walls. The effect of differences in kinetic rates is also expected to decrease 
as the flow field becomes more reactive, i.e., at higher stagnation pressures. 

A high density, small nose radius comparison is shown in Figure 5-3. Under 
these conditions the effect of the two different sets of reaction rates is 
negligible . 





CONDITIONS AND RESULTS FOR COMPARISON CASES 



for N and 0, zero for NO 



























































FIGURE. 5-*2- EFFECT of HOMCq£N£OU 6 FIGURE S-3 Sr-FSCT OF WOMOiSENtEOUS tdlWETlCS 

KlNJETlCS 










5.2 


effect of numbfr of nodal point 


The BLIMP computational procedures uses spline fits of the primary 
variables across the boundary layer and therefore require a lesser number of 
nodes than linearized finite difference methods. Appropriate nodal point dis- 
tributions for equilibrium chemistry are recommended in Reference 1 and it 
appears logical to determine if similar distributions are valid for nonequilib- 
rium conditions. A typical nodal distribution has 7 nodes for laminar boundary 
layers and 13 nodes for turbulent boundary layers. These distributions are 
not sacred, however, they and similar distributions provide good starting 
points for comparing the effect of nodal distribution on predicted boundary layer 
variables. In this and most of the subsequent comparisons, the mass fractions of 
N, 0, and NO will be shown and used as criteria for determining the magnitude 
of any discrepancies. 

The stagnation point specie distributions for a catalytic wall and a 
near-noncatalytic wall are shown in Figures 5-4 and 5-5, respectively. The 
solid lines represent a 13 node distribution and the indicated "data" points are 
from a 7 node distribution. Both distributions are shown in Figure 5-1. Also 
shown in Figure 5-4 is the fqui librium distribution. The strong resemblance 
between the two profiles in Figure 5-4 is a result of the high catalycity which 
reduc ;s atomic species to zero at the wall. Between the 7 and 13 node distribu- 
tion there is a small but noticable discrepancy in the NO concentrations for a 
catalytic wall. However, this is accented somewhat by the fact that the NO 
scale is magnified by an order of magnitude. All differences between 7 and 13 
node distributions are considered small and result in no significant difference 
in predicted heat transfer ratet as shown in Table 5-1. 

5.3 EFFECT OF SURFACE CATALYCITY 

The effect of catalycity on the specie distribution within the boundary 
layer are shown in Figures 5-6 to 5-8. It was assumed that the surface catalytic 
reactions were given by 
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The catalycity for all three reactions were assumed equal although it will be 
subsequently shown in Section 5.6 that the apparent catalycity of typical shuttle 
materials is greater for 0 recombination than for N recombination . Finally, it 
was assumed that all reactions were first order so that the catalycity can be 
defined as 


K (A' = p D [SP (A) ] 
w w w 3y 

where A represents 0, N or NO and SP(A) represents, the mass fraction of A. 

It can be seen from Figures 5-6, 5-7 and 5-8 that the extent to which 
the reactions at the wall influence the boundary distribution is much greater 
for N than for o. This is a result of the strong influence of the nitric oxide 
shuffle reactions which has the net effect of wanting to deplete N in favor of 
0. r r , 0 atom concentrations even for highly catalytic walls are virtually 
frozen for about 3/4 of the boundary layer. This effect will be shown in 
Section 5-4 to decrease (i.e., a smaller fraction of the boundary layer being 
frozen) as the pressure decreases and vice versa. 

The ratio q/q is shown in Figure 5-9 as a function of the catalycity 
cat 

for the above cases which correspond to a velocity and altitude of about 
21,600 ft/sec rnd 225, 0C0 ft, respectively. This is compared with the results 
of Reference 40 which used a si aller nose radius (1.95 ft as opposed to 2.325) and 
a lower vail temperature (7C0°K as opposed to 1800°K) . It can be seen from 
Figure 5-. that a significant difference exists between the two predictions. 

The effect of the difference in nose radius is expected to be small. An increase 
in wall temperature would shift tl.e curves to the right slightly, however, most of 
the difference is probably due to the binary gas behavior used in Reference 40. 
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It is clear from Figure 5-7 that NO distributions are strongly dependent on the 
wall catalycity and the resulting interaction with the shuffle reaction will 
doubtless have significant effects on the wall values of 0 and N concentrations 
(thereby affecting the heating rates) . It is, therefore, necessary to exercise 
some caution in the specification of wall catalycities from experimental data 
since homogeneous reactions, even in near frozen flows, are intimately coupled 
to the wall reactions. 

5.4 EFFECT OF PRESSURE 

The effect of stagnation pressure on specie distributions is shown in 
Figures 5-10 and 5-11 for a low catalycity wall (144 cm/sec) and a high catalycity 
wall (3600 cm/sec) . Since the boundary layer thickness is highly dependent upon 
the pressure the normal cordinate y was normalized with respect to a reference 
value representing the value of y at a fixed value of u/ue (in this case u/ue - 
0.8) . The degree of dissociation at the boundary layer edge for equilibrium is 
also dependent on the pressure, however, no normalization was made in plotting 
the graphs. Although there are differences in the distributions at different 
pressures, these differences are not dramatic a= long as the wall catalycity is 
fixed. 

A comparison of the heat transfer rates as a function of pressure are 
shown in Figure 5-12. Also shown are the slopes for q " /"jp” . As might be 
expected, catalytic wall heat transfer can be scaled by this classic relationship 
with moderate accuracy, however, low catalycity walls can not be scaled in this 
form. 


5.5 COMPARISON Ov CURRENT PREDICTIONS WITH REFERENCE 61 

The relative accuracy of predictions from the current code were compared 
with those from the computer code described in Reference 63 . This latter code 
was made available to Aerotherm by F. G. Blottner so that ou^-mt from both codes 
could be compared directly, A case from Reference 61 was used for comparison, 
namely, a spherical nose with a rad’ us of n .0833 feet at stagnation conditions 
of 6.026 atmospheres and 7886 Btu/lbm. The rate constants used in both codes 
were those shown in column 1 or Table 5-2 a.-.u edge conditions were those given 
in Reference 61. Specie distributions are shewn xn Figure 5-13 for a wall which 
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is catalytic to N and 0 recombination, but noncatalytic to NO reactions. Similar 
distributions are shown in Figure 5-14 for a noncatalytic wall. For catalytic 
walls, the ditferences are relatively small, however, a significant difference, 
especially in the concentration of 0, is noted for the noncatalytic case. 

The probable reasons for these differences have not been resolved but would ap- 
pear to be either in the computational procedure or the thermodynamic models. 

5.6 ANALYSIS OF ARC JF.T DATA 

Plasma arc jets have been used extensively for the simulation of reentry 
heating environments and some recent data was reported in Reference 62 for 
candidate shuttle materials, namely, oxidation inhibited carbon-carbon composites 
(LTV) , coated silica refractory insulation (LI-900) and coated coluxnbiun. The 
surface coating on the carbon-carbon is believed to be primarily silicon carbide. 
Composition of the other coatings is given in Table 5-3. All coatings have a 
significant amount of silica which at ligh temperatures is expected to form a 
glassy surface with a low to moderate catalycity. The referenced data was 
obtained for dissociated air over a range of enthalpies between 3000 and 20,000 
Btu/lbm and at a pressure of about 0.01 atm with flat faced cylinder models. 

At the low enthalpy conditions, only oxygen would be dissociated under equilibrium 
conditions whereas at the higher enthalpies both oxygen and nitrogen would be 
dissociated. The maximum enthalpy of 20,000 Btu/lbm is significantly greater than 
would be experienced by an earth orbit reentry vehicle and would result in some 
ionization which is not considered in the current analysis. The data was not 
obtained with the intent of calculating catalytic efficiencies so that there is 
no data for any one of the materials which span the whole enthalpy range. However, 
since all coatings contain silica and nave glassy properties, it was assumed that 
their catalytic behavior would be similar. Thus, the data was considered as a 
complete set and the 3LIMP/KINET code used to estimate the average catalytic 
efficiencies for both O and N recombination in dissociated air. 

The calculations were performed on the assumption that the edge of the bound- 
ary layer at the stagnation point was in chemical equilibrium. Although the chemi- 
cal state of the arc jet flow was not measured, the stagnation point would be at or 
close to equilibrium since a high pressure dissociated equilibrium state is at- 
tained in the plenum and the rapid expansion in the nozzle retards chemical recom- 
bination. Thus only a small change in chemical composition would ba required, as 
the flow transverses the bow shock, to reach a stagnation point equilibrium state. 
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TABLE 5-3 

ELEMENTAL COMPOSITION OF COATINGS, % BY WEIGHT 



R5 12E 
Coated 
Columbium 
752 

VH 109 
Coated 
Columbium 

LTV 

Coated 

Carbon- 

Carbon 

LI 900 

Si 

38 

70 



Nb 

57 




W 

4 




Fe 

.2 

5 



Cr 

.1 

5 



Hf 


20 

0) 

(2) 


(1) Proprietary but believed to be mostly SiC 

(2) Silicon carbide, borosilicate glass. Ref. 22 
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Stagnation point solutions were then obtained over the prescribed enthalpy 
range for various values of the surface catalycity for O and N recombination. For 
all but one set of computations, the Wall was assumed to have no effect on NO 
chemistry. The assignment of a wall that is no-'catalytic to NO reactions was 
simply an expedient to reduce wall catalytic effects to two parameters, namely, 
the catalycities of N and 0. It is clear that a significant number of coupled 
chemical reactions occur simultaneously on the wall and in a zone near the wall. 
Further it is known that the gas phase NO shufiJe reactions are very rapid so 
that slow generation or depletion rates of NO at the wall would be rapidly compen- 
sated for by homogeneous gas phase reactions. Even s , since it is not possible 
to determine the r.tes for all possible simultaneous surface reactions, the data, 
from which multi-component catalycities are deduced, should be obtained at pressure 
levels which are representative of the flight environment. 

xhe wall reaction rate parameters (FKF) are shown in Table 5-4 for each set 
of calculations and the. ratio of q/q is shown in Figure 5-15. Also shown are 

C3t 

the data for the three materials as determined from the equilibrium radiation sur- 
face temperature. At low enthalpies where only oxygen is dissociated, the curves 
all come together at a value of FKF (0) equal to 0.02. At higher enthalpies where 
both oxygen and nitrogen is dissociated, extensions of the curve for FKF (N) equal 
to 0.006 appear to best represent the data, at least for enthalpies less than 
14,000 Btu/lbm. At higher enthalpies the ratio q/q either remains constant or 

C3t 

increases slightly. This may be due to ionization effects or a shift in composi- 
tion of surface adsorbed atoms from 0 to N, thereby increasing the catalytic ef- 
ficiency for N recombination. At the higher enthalpies, due to the fact that ioniza- 
tion was not included and complete dissociation was achieved, the amount of energy 
stored as chemical energy no longer increased with the total enthalpy of the rlow. 
Thus, fortuitously, the predicted ratios of q/ c 3 ca t also increased at the higher 
enthalpies. 


The values of FKF <\re related to K w (as shown in the Appendix) by the 
relationship 


T 

FKF x w 

K W “ 0.761 


ft/sec 


when the surface reaction is considered as 


+ ~ A 
2 2 
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TABLE 5-4 


SURFACE REACTION RATE CONSTANTS 


URVE 

FKF(O) 

FKF(N) 

FKF(NO) 

1 

O.C2 

0.006 

NC 

2 

0.02 ‘ 

0.012 

NC 

3 

0.02 

NC 

NC 

4 

0.02 

0.02 

0.02 

5 

0.02 

0.0002 

NC 

6 

0.02 

0.002 

NC 

7 

0.002 

0.002 

NC 

8 

0.02 

0.008 

NC 
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SECTION 6 


APPROXIMATION OF NONEQUILIBRIUM CONDITIONS AT 
EDGE OF BOUNDARY LAYER 

The development of a nonequilibrium computer code for calculating the 
inviscid flow field ab^ut a three-dimensional shuttle body is a formidable task 
being pursued by several investigators. Until these codes become operational, 
however, some approximate methods are required to obtain realistic boundary 
layer calculations. In this section, an approximate procedure is preserted for 
specifying the edge boundary conditions. The computations are performed in three 
steps. The first step is an approximation of the pressure history along various 
streamlines in the inviscid flow field. The second step uses the pressure his- 
tory and initial conditions of the streamline to calculate the thermochemical 
state, including entropy, along the streamline. The final step is to match the 
boundary layer mass flow with the mass flow represented by the streamline. 

Clearly planar or axisymmctric assumptions are necessary in order that mass flows 
be defined. 

Although in the current context, this procedure is applied to shuttle 
flow fields, it is not necessarily limited to this application. For instance, 
the chemical relaxation code (step 2) can be used for the computation of chemi- 
cal reactions on any prescribed pressure streamline. 

G.l CORRELATION OF INVISCID SHOCK LAYER FLOW FIELD 

The inviscid flow field on the windward pitch plane of a space shuttle 
was calculated over a range of flight conditions representative of a typical 
entry trajectory. In particular, the initial conditions immediately behind the 
detached shock wave and the spatial pressure history of several streamlines 
were calculated. The correlations presented here were the data base for these 
calculations . 

The crucial features known to dominate the flow field property distribu- 
tion are the shock surface and the body which supports it in a hypersonic ang^e 
of attack flow. For this reasc , curve fits were prepared for second-order (cur- 
vature variation) smooth shoe): and body surfaces and used as defining boundaries 
in the correlation functions. 


\ 

J 




l 



In addition it is convenient to characterize the equilibrium hypersonic 
shock layer equation of state properties using an analytic representation as 
opposed to more cumbersome table look-up. These state functions are also pre- 
sented he.:e. 

The data correlations are supplied in graph and tabular form and where 
available, in functional form for subsequent use. Included are pressure distri- 
butions, enthalpy ratios at the shock entrance position (transition ratios) for 
each of the selected streamlines, streamline traces (distance along the stream- 
line as a function of body axial coordinate station) and quantitative comparisons 
with available measurements or exact inviscid calculations. 

6 . 2 GEOMETRY 

Figure 6-1 shows the body oriented coordinate system adopted for the code 
solutions as well as the correlations. Here the windward generator axis is 
rotated about its origin at positive angle of attack, a. The axial distance is 
x measured from the origin It = 0 at nose. The basic coordinates are the distance 
along the body measured from the forward stagnation point, £ Q , or along selected 
streamlines, measured where they enter the shock, £^ = 0, for example, which is 
the entry position for the ith streamlire. The streamlines selected for the 
correlations are traced through the shock layer a distance £ which corresponds 
to a body station distance (axi ’ dimension) x, given in the graph. Figure 6-2. 

All length measures in the correlations are ratioed to the nose radius, 

R q . r is the dimensionless radial dimension of a point in the shock layer meas- 
ured from the body generator axis while r w (x) is the radial surface dimension 
at station x, and r^ (x) is the radial shock dimension at x. y is the coordinate 
measure of a point in the shock layer measured along a ..ormal from the body and 

has a value A (the shock standoff) at the shock. 6, and 6 refer to shock and 

o w 

body inclination angles, respectively. Both are measured as clockwise angular 
deviations from the normal to the body axis as shown. The effective shock angle 
at positive angle of attack, a, is simply 

AH (0 6 - a) • (150) 

The components of velocity, u, v are oriented positive along streamlines 
£ and outwari normal to them, respectively. Velocities and components appear as 
ratios to the free stream velocity, U a . Shock layer density p” is ratioed to 
the free stream density p M . Shock layer pressure p is ratioed to twice the free 
stream dynamic pressure, although in the correlations, it always appears as ra- 
tioed to the stagnation point pressure P To » or the local streamline total pres- 
sure, P T< 5 > The sta.ic enthalpy h is ratioed to the square of the free stream 
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velocity, although in the correlations it is normally ratioed to the total 
enthalpy, H^. 

Figure 6-3 shows details of the shock layer stand-off geometry. Once 
the body surface variables r (x) , 9 (x) , < (x) are computed from the curve fit 

WWW 

expressions, the necessary shock values may then be readily computed from the 
wall variables and one additional correlatable parameter, the shock stand-off, 
6(x) The geometry of Figure 6-3 is characterized by a simple set of defining 
shock wave-surface relations: 



( 151 ) 

(152) 

(153) 




+ dr. 



(154) 

(155) 


The mass flow calculations which are used to obtain the y displacements 
of the streamlines in the shock layer, and to establish the position at which 
a streamline enters the boundary layer are based on some additional geometric 
considerations . 

Figure 6-4 shows the geometry illustrating the mass flow balance used 
for the inviscid streamline displacements. Adoption of a generator axis nose 
intercept origin for the axial displacement (x = 0) necessitates accounting for 
both the radial offset of the generator axis (A in the figure) and the normal 
dimension, S q , which is a measure of the shift of the forward stagnation point 
at angle of attack. The mass flow intercepted by the shock over an annular 
section of radius, S q , is spilled over the leeward side of the body at angle 
of attack. This is accounted for by subtracting this leeward mass flow when 
computing the dimensionless mass flow for a given streamline. 


6-5 



I 



FIGUEE. <^-3 -6 HOCK 5T AMD-OFF GEOMETRY 



FKqurae. g-4 


MASS FLO\M 6£OMeT(?t 


6 -G 


The geometry leads to the following set of relations for dimensionless 

A 

mass flow associated with the ith streamline. In these equations Z is the 
rotated angle of attack annular radius while y is the local flow angle. 

Z. = (r. + A.) cos Qt (156) 

1 01 

A^ = (x^tan a) (157) 

in. = (Z 2 - S 2 ) = n[{ (r . (x) + A.) cos y} 2 - { (r (x) + A) cos y} 2 ]pu cos y (158) 

1 1 O 1 1 w . 

1 

6.3 CASES 

A selected number of cases from a typical shuttle trajectory (see Section 7) 
were covered in the correlations and subsequent calculations. Flight conditions 
for these cases, which span a Mach number range from 9.61 to 30.2 and an angle-of- 
attack range from 30 to 34°, are shown in Table 6-1. Although the correlations were 
obtained for a limited angle-of-attack range they are expected to be valid over a 
much wider range. 


TABLE 6-1 

NAR BASELINE TRAJECTORY 


Case 

Number 

U 

00 

(kfps) 

Altitude 

(kft) 

M 

00 

O 

a 

t 

(sec) 

1 

25.6 

300 

30.2 

34 

250 

2 

25.4 

250 

28 

34 

400 

3 

23.9 

238 

25.1 

33.3 

600 

4 

21.8 

224 

22 

32.3 

800 

5 

18.7 

201 

17.9 

31.1 

1000 

6 

14.6 

181 

13.6 

30.2 

1200 

7 

10.4 

164 

9.61 

30 

1400 
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6.4 


EQUATION OF STATE 


The state relations for the equilibrium air shock layer calculations were 
correlated from the normal shock parameters given in the Cornell Aero Lab tables 
prepared by Witliff and Curtiss (Ref. 63). These are based on the 1959 ARDC Model 
Atmosphere . 

The equation of state model is a particularly simple form yielding sur- 
prising accuracy in hypersonic equilibrium blunt body shock layer calculations 
of many investigators. In particular we use the form suggested by Swigart (Ref- 
erence 64) in his 1963 theoretical investigation of the blunt body angle of attack 
problem in equilibrium air 


h = 



A 

o 


h = h/u 2 P = p/p u 1 

00 r CO CO 

= h/u 2 o' = p/p 

0 0 00 ^oo 


(159) 


where y is an effective caloric ratio on a particular adiabat (shock entry point) 
and A^ is the energy deviation from polytropic due to equilibrium gas phase 
chemistry. 

The current range of calculations correlate within about 5 percent maximum 
deviation from the shock table values to the following forms 

A q = .098 * (M*)~ ,77 , 12£M*£31 (160) 

M* < 12 

OO 

Y =1.5* (M*)"* 085 , 9.5£M*£31 (161) 
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by 


The effective shock Mach number, M* for the angle of attack cases is given 


M* = M cos A 
00 00 


(162) 


These shock state correlations are altitude invariant over the range of 
altitudes treated in the 7 cases (Table 6-1) . The shock Mach number range of 
applicability is recorded with Equations (160) and (161) . 


6.5 BODY SURFACE AND SHOCK SURFACE REPRESENTATIONS 

To insure smoothness in both local body surface slope and curvature, the 
body surface was represented by a connected sequence of five second order poly- 
nomials. The representation is initiated at the generator axis intercept at the 
shuttle nose (x = 0) . Avoidance of higher order polynomial representations or 
alternate functionals insured timely smooth surface profile and derivative gene- 
ration without embarking on a lengthy process of tedious analytical smoothing. 
The body generator equation in the windward pitch plane is 

r" 2 = a + bx + ex' 2 (163) 

w 


TABLE 6-2 
EODY COEFFICIENTS 


Section 

X 

a 

b 

c 

I 

[0, 0.525] 

0 . 

1.92059 

-0.84868 

II 

]0.o25, 2.60] 

0.2550 

0.94915 

0.07651 

III 

]2.60, 15. 2[ 

-0.40705 

1.45843 

-0.02143 

IV 

[15.2, 46. 5[ 

3.11583 

0.99502 

-0.00619 

/ 

[46.5, 80[ 

r (46. 5) 2 
w 

0 . 

0 . 

. 


(cylinder) 
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The symbols [a, ,b[ indicate whether or not the range limits a,b are 

included or excluded from the interval, respectively. No attempt was made to fit 
the boat tail region. 

The body and shock curves were based on extensive wind tunnel shadowgraph 
measurements of shuttle profiles at angle of attack reported in previous Aero- 
therm analysis (Reference 65) . These data were supplemented by experimental 
shadowgraph traces reported by Marvin et al. (Reference 66) and numerically exact 
inviscid angle of attack model shuttle calculations presented by Rakich and 
Kutler (Reference 67) . 

A particularly significant factor in the development of the shock shape 
correlation for angle of attack shuttle flow fields was the observation that 
the geometrical relationships between body and shock surface were invariant de- 
spite changes to the angle of attack. The angle of attack effects on shock 
transition variables between about 25° and 50° positive pitch were recoverable 
by rotation of an assumed fixed shock surface-body surface unit. The shock 
stand- wff adjustment was correlated by accounting for Mach number change with 
density jump across the rotated shock. 

A set of relations for predicting the local shock stand-off at a specified 
body station, together with the shock shape equations, (151) through (155) , are 
sufficient to generate the shock surface of a given Mach number, altitude and 
angle of attack of interest. The stand-off correlations are as follows. 

In the "afterbody" region (x >_ 1.) the ratio of the "baseline" (M^ - 9.61, 
a = 30° , altitude = 164 kft) to its value at the generator axis is computed in 
three segments: 

(6/6 ) = f(x) =0.139 + 0.891 x 1. < x < 15. 

Bo — 

(6" /6 ) = f(x) = 7.7 + 0.387 x 15 < x < 46.5 (164) 

BO ' 

(6/6* ) = f(x) =-10.24 + 0.773 x 46.5 < x 
B O 
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The baseline "deltas" may then be corrected for the particular Mac^ 
ber altitude and angle of attack 





The effective shock angle. A, was introduced previously. 

In the foregoing Equation (165) the dimensionless normal shock wave 

stand-offs, 6 and 5 (for the particular and baseline cases) are computed from 
O BO 

the hypersonic similarity spherical shock expansion written in terms of the in- 
verse of the normal shock wave density ratio, e 


4 > 

T » ■ -*T- 

<j>£ = e - e 3 ^ 2 + 3e 2 + a(e 3 ) 


(166) 


In the forebody region the shock stand-off variatic. from its value at 
the generator axis, o q , to its value at the shifted stagnation point, , on the 
windward side is correlated by combining che results of Kaatta’-i (Reference 68) 
and Swigart (Reference 64) . First compute the radial position of the forward 
stagnation point 


(r )__ = 1.34 sin a - 0.682 sin 2,2 a 
w oP 


(167) 


Next compute the variation in stand-off from r = 0 to (r ) 

W W SP 


H = 6 - (r 

o 


“ ° u <a7 


0 < r < (r )__ 
w — w SP 


(168) 


o w 
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where 


1 j 

— “ = 0.94 sin a - 0.046 
dr 

o w 

Over the short distance between forebody stagnation point and the after- 
body range correlated by Equation (164) , a linear growth of stand-off correlates 
the data 

6 - « SP + («<1> - & et ) ( ) (169) 

\1 - x gp / 

Two additional dimensionless metrics, the shock wave normal, N(x) and the 
shock wave tangent metric, C, are used in the shock relations developed subsequently. 
Both depend on the derivative of the shock wave trace 3x/8'r^ , as does the local 
shock slope. This derivative is obtained as follows. Combine Equations (if j 

(152) obtaining 


x = x + 6 - cot 0 (r x - r ) 

o w 6 w 


Differentiate and rearrange, obtaining 


d0 


/ dr. 


X* = ( — )= x' + (r,, - r )cso 2 0 [ — — 1+ cot-0 ( — — - 1 


3r, 


l r. r- i 

o w 


w \ — 
.dr. 


W \ 


(170) 


dr 


w 


Differentiate Equation (151) and neglecting terms 0(e) or smaller we obtain 


dr 


w 


= 1 + 6(e) 


(171) 


dr, 


and d0 /dr can be obtained directly 
w w 
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<30 /dO \ / dC \ ,r ( 

#' ir F 

w \ / \ w / V 


1 + (X 1 ) 


• 1 2 


(172) 


Introducing (171) and (172) into (170) the expression used for the re- 
quired derivative is obtained 


x' = x' 


+ K (r - r ) \Jl + 


(x') 2 / sin 2 6 


w w 


w 


(173) 


For 0 > 0 

w 

and x' =0, for 0 =0 

w 

The body slope, x' and the body curvature are obtained directly from 
differentiating Equation (163) . 

A comparison of the accuracy of the shock and body correlation vs measured 
values is presented in the following graphs. 

Figures 6-5, 6-6, and 6-7 show comparisons of computed 'S measured 

* 

(shadowgraph) traces of the shock slope supported by the shuttle model bodies 
(Refc-ncer PI and P3) . The more sensitive body and shock slope results are com- 
pared in Figure 6-7. Figures 6-5 and 6-6 show the overall results and details 
restricted to the forebody region respectively. Profile 2, is the selected body 
shape function generated by the coefficients in Table 6-1. 

6.6 COMPUTATION OF THE SHOCK TRANSITION VARIABLES 

In the generalized body-oriented coordinate systems for the shock layer 
previously introduced, a formal development of the jump relations is obtained. The 
metric direction normal and tangential to the shock surface are readily determined 
in terms of the foregoing body and s''ock surface rel*. ti is and their derivatives. 

With this information, the remaining shock transition variables are readily 
computed for shock intercepts in terms of arbitrary body stations x, selected for 
th 2 analysis . 

Formulation of the shock relations are an extension of the previously re- 
ported NASA supported work of References 69 and 70 with modifications suggested by the 
formulation of Webb et al., (Reference 71). 
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Reference is made to the principal axis vector direction cosines shewn in * 
Figure 6-8, the geometry introduced in Figures 6-1 and 6-2 and the non-dimensional ^ 
scheme used for the flow quantities, introduced in Section 2. 

H 



Figure 6-8 Direction Cosines 


We choose e^ = - e^, where $ is the aximathal angle and the direction 
metric is outward. 

The direction cosines are written 

(body normal) = e x COS °b ~ e y Sln r 


(174) 


(body tangent) 
(binormal) 


t, = e„ sin 9. + e„ cos 9 
b X b Y b 


b = - e = - e. 
b Z cj> 


Multiplying by cos 0^ and t^ by sin 0^ and adding, yields fox; e x 


= r, cos 9,, + t^ sin 0. 


'X b 


(175) 


Multiplying n^ by -sin 0^ amd t^ by cos 0^ and adding , yields for e^ 
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e„ = - r, sin 6. + t. cos 0. 
Y b b b r 


(176) 


While, quite simply 


-y -y 

* ' b b 


(177) 


The normal and tangent surface metrics of a shock surfece X(rg,4>) arbi- 
trarily oriented in our body geometry for positive angle of attack may then be 
written 


”*( 


.. 2 W2 

1 + X' 


C = cos 


a + sin a ^ x* ^ 


(178) 


The derivative of the shock trace is obtained using Equation (173) , 
the <j> derivatives are neglected close to the plane of symmetry on the windward 
side. 

The free stream velocity can next be decomposed into components exterior 
and tangent (t^) or normal (n^) to the shock surface at angle of attack, under 
the pitch plane symmetry simplification. See, for example. Reference 71. 

The components are 


(U ) 

« n c 




(179) 


% = 
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To obtain the relations used first write the velocity component directions 
(179) and (180) in terms of the body direction cosines, equations (175), (176), 
ar i (177) , obtaining 



Now ' .e free stream vector velocity may be expressed compactly in terms 
of the developed normal and tangential metric coefficients. Equations (181) and 
(182) . 


Uoo = ( u -|)^6 + { u ^ 1 ~ c 2 / h 2 ) : 


(183) 


' e Rankine-Hugoniot relations for an oblique shock wave are summarized 


as 


Ur. = £u 
6 '*■ 


h 6 * h ~ + f (1 ' e 5 


V. = V 

<5 CO 

Pr = p + p u 2 (1 - e) 

*0 c oo ^OO 0O 



(184) 
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For computation we combine the foregoing metric expression for the free 
stream velocity, Equation (183), the Rankine-Hugonict relations, Equation (184) 
and the equation of state. Equation (159) , obtaining the relations in order of 
their calculation. For any selected station x, following calculation of all of 
the body surface ar hock surface variables compute 




(188) 


(189) 


and 


G 


2 


= _x_ 

Y-l 



In addition we compute the total pressure at the local streamline inter- 
sect with the shock front 




(S 6 + * 6 > 


(190) 


and the total enthalpy 


H. 


t 6 - h 6 + I + ^ 


(191) 
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Using these to form the streamline to streamline correlation ratios 


P 




Equations (185) through (191) , the shock surface and body surface relations 
and a streamtube integration of the streamwise momentum were calculated for the 7 
trajectory points using a desk calculator and the correlation relations for 9 
streamlines. 


6.7 STREAMTUBE CHEMICAL RELAXATION 

The streamtube code was developed to interface with the BLIMP/KINET code 
and the prescribed pressure distribution presented in previous sections. By 
matching mass flow rates it is thus possible to approximate both nonequilibrium 
chemistry and variable entropy at the edge of the boundary layer. Although any 
available reacting streamtube code could be used, the code described below was 
written so that the calculation of thermochemical data would be identical to that 
used in BLIMP/KINET. Thus there would be no incompatibilities m the chemistry. The 
computational procedure is implicit and is therefore numerically stable even for 
large step sizes. 

The equations to be solved are the combined energy and momentum equation 


dh _ 1 dp 
ds p ds 


(192) 


and a set of n-1 specie equations 


d_ 

ds 


(SP i ) 


pu r i 


(193) 


subject to the constraints that 
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- h = 


u 


2 


2 



( 194 ) 


n 



i=] 


SP. 

r 


« 1 


Let m denote the spatial station along the streaintube and £ be the iteration at 
station m. Then a Newton-P-aphson procedure was used to solve the above equa- 
tions. Assume that the solution is known at m-1 (which could for instance be 
the initial conditions) then the independent variables h and SP^ at station m 
for the &-th iteration is related to the known m-1 station and S.-1 iteration 
by the equations 


s 


H 


I 



( 196 ) 
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where 


Ah = h o - h „ . 
m,£, m,x,-l 


Asp. = (sp.) _ - (sp ) . . 

i x m,x, x m,x.-l 


As = s - s , 
m m-1 


( 197 ) 


HHH 


m-1 


_ L + As / 1 §£ \ ^ 

- Vl^pdsJ^j 


As { 




Equations (195) and (196) represent n equations with an equal number of unknowns, 
namely Ah, ASP^ ASP 2 * . .ASP^ . Implicit unknowns are 


p = p(SP i ,h,p) 


T = T(SP i ,h,p) 


(198) 


<P i = <J> i (T,SP i ,p) 


u = u (h) 

Equations (195) and (196) are solved iteratively with a step size constraint 
based simply on the number of iterations required to obtain a converged solution. 


« 
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SECTION 7 


CALCULATIONS FOR REPRESENTATIVE SHUTTLE VEHICLE 


For these calculations nonequilibrium chemistry in both the inviscid and 
viscous regions were considered. Pitch plane calculations were made for the 
windward side of the Rockwell International orbiter (Figure 7-1) at flight con- 
ditions representative of their 2007 baseline trajectory (Figure 7-2) . The pitch 
plane outline of the vehicle at angle of attack was assumed to be the generator 
of an axisymmetric body with a shock shape determined from correlations of wind 
tunnel shadowgraphs (Reference 65) . 

A simple-in-theory but involved-in-practice procedure was used to estimate 
the nonequi] ibrium and nonisentropic conditions at the edge of the boundary layer. 
The following steps were required for each trajectory point. 

1. Calculate body pressure distribution 

2. For a prescribed shock shape, estimate the pressure distribution along 
1 various streamlines in the inviscid flow. 

3. For each streamline, starting with a frozen oblique shock state, 
chemically relax along the prescribed pressure-distance history to 
yield the thermochemical state along each streamline. 

4. Estimate the boundary layer mass flow using equilibrium assumptions 
and compare with the mass flow represented by each streamline to deter- 
mine where that streamline should enter the boundary layer. Then from 
3 the thermochemical state at the edge of the boundary layer is known. 

The current inviscid flow calculations are not applicable in the stagnation 
region except in an integrated sense. In fact no streamline approach would be 
valid on the stagnation line for nonequilibrium conditions because of the zero 
velocity limit. Thus, the nonequilibrium solution in the stagnation region was 
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calculated using a viscous shock layer option and these solutions were then matched 
with the solution from (4) above to use as the edge conditions for a boundary layer 
calculation. 

7.1 SURFACE PRESSURE DISTRIBUTIONS AND APPROXIMATE BOUNDARY LAYER MASS FLOW RATE 

Surface pressure distributions for the cases given in Table 6-1 were 
computed using the modified correlations presented in Reference 71a and are shown 
in Figure 7-3. Similar distributions could be calculated via a smooth transition 
from modified Newtonian values in the vicinity o the nosa to local tangent cone 
values downstream as suggested in Reference 65. With these pressure distriuutions , 
boundary layer calculations* were used to determine the bour lary layer characteristics 
on an equivalent axisymmetric body. Since axial symmetry was assumed, the mass 
flow rate at any given body station is defined and can be related to a stream tube 
mass flow. Typical mass flow rates for 4 cases are shown in Figure 7-4. Of the 
cases given in Table 6-1, only cases 1, 2, 4 and 7 (t = 250, 400, 800, 3400 seconds) 
were considered for complete heating analysis. 
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Since chemical reactions along an inviscid stagnation streamline are not 
well defined, the BLIMP/KINET viscous shock layer option was used to determine the 
thermochemical state in the shock layer. A nor.catalytic well was assumed since 
primary interest was in the state at some point which would represent the edge of 
the boundary layer. It was also assumed that the shock wave was chemically frozen. 
A typical set of distributions for t = 800 seconds is shown in Figure 7-5 through 
7-7 for three body stations. The rapid dissociation immediately downstream of 
the shock wave is evident however the remainder of the shock layer is not very 
reactive. Similar shock layer solutions were obtained for t = 250, 400, and 1400 
seconds for estimating Lo'indary conditions as described in Section 7-4. 


7.3 STREAMLINE PRESSURE AND ENTHALPY CALCULATIONS 

The enthalpy, pressure, density, and velocity components were computed 
for 9 selected streamline entrances at the shock front including the body streamline. 
These streamline shock jump conditions and a determination of the pressure expansion 

♦Integral boundary layer calculations wc,.e used here, however any reliable procedure 
including the BLIMP code could be used. 
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downstream on individual stream tubes provided the initial and boundary conditions 
for the subsequent flow field chemistry and reacting boundary layer solutions. 

The jump conditions are computed at arbitrarily selected positions on the 
computed shock surface from the shock relations developed in the previous sect; ->r.- 
In the presented analysis, the ratios P_/P «, H /H , p , u , and v were found to 

4 1 4 1 4 4 4 

be useful parameters for application to the subsequent code solutions. 

From the shock wave surface relations, equations 151 through 158 and Equation 
(173), the parameters 6, , 0^, X^, X^"*, r^ are obtained. Next the initial sweep 

or arc of the streamline as it departs from the shock correlates simply as the 
sector of a characteristic circumscribed semicircle in the forebody (and tends to 
zero in the afterbody) 


, _ TTRpS/ R 

5 2 


The normal displacement. 


(199) 


Y * 



of any streamline is next computed for all streamlines which entered the shock up- 
stream of the shock station being computed. This is accomplished using Equation 
(158) to determine the dimensionless mass flow, m for all streamlines. The pro- 
cedure is, as follows. Consider a streamline entering the shock at the jth station, 
the upstream entering streamlines (J-l, J-2 , ...) have normal displacements, measured 
at the jth station, given by 


III = / 


m j-l /m j 






(200) 


The streamline path lengths are then computed, 

J 


tt6_ m 
(S ) I = -Jill i 

U 


K 


^ ' Vl 1 ’ + <*K ' 

— u — N 


(201) 
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A plot of the calculated streamline path lengths as a function oi oody sta- 
tion axial coordinate X appears in a previous section (Figure 6-2) . 

To compute the pressure distribution along the streamlines correlative 
procedures were developed. Observations derived from wind tunnel data and exact 
numerical inviscid flow field calculations were used as a basis for the pxroi. lure 
and relations developed. 

A regularity in the body surface pressure distribution for both adiabatic 
and nor-adiabatic shock layers supported by long, analytically smooth, blunt- 
nosed body surface profiles at zero angle of attack has been exploited in many 
inviscid flow analyses. In particular, under a NASA sponsored study Buckingham 
and Hoshizaki (Reference 72) showed that for affinely related power law bodies 
with slenderness ratios (L/D > 3) series correlations exist for both pressure 
and convective heat transfer that are independent of both Mach number and 
altitude. These correlations are verifiable for a restricted but useful range 
of hypersonic reentry trajectories. The profiles, which posses spherical, 
oblate and prolate ogival, or paraboloid noses generate correlatable pressure 
distributions. Success of the correlations depended, in part, on treating the 
pressure distribution as a ratio of the local pressure to the total pressure on 
the streamline, P/P^. 

From studies of exact numerical inviscid flows on shuttle vehicles at angle 
of attack (Rakich and Kutler, Reference 67) experimental wind tunnel shuttle data 
(Marvin et al Reference 66) and the extensive data tabulated by Bartlett, Morse, 
and Tong at Aerotherm (Reference 65) , the constancy of the ratio P/P T on each 
streamline along normal y at a particular body station may be noted. Some results 
of previous theoretical studies, using exact numerical method of characteristics 
and finite difference solutions help to substantiate this observation. For in- 
stance, Figures 7-8 and 7-9 are for axisymmetric flow past an ogive (L/D =3.5) 
and sphere -cone -cylinder (L/D = 12) method of characteristics with equilibrium 
air (Reference 73) . Similar results have been obtained by reducing the profile 
data of Rakich and Kutler (Reference 67) and experimental data presented by 
Marvin (Reference 66) . Significant deviation from constant P/P T2 along a normal 
was restricted to the forebody region of this representative sample of both 
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shuttle and shuttle-like slender bodies of revolution at various angles of 
attack. To anticipate these results examine the X and Y inviscid momentum 
equations. For the adopted coordinate system, in dimensionless terms, these 
are: 


— Su ,, — . — 3u 

P u »r+ (l + tcy) pv ~ = 


3p 

3C 


a, 

'M) 


'VO 


( 202 ) 


— g v — 3v <pu 2 

pu ac + pv ■ (i + Ky> “ 

% % \ 


9p 

3y 


(203) 


Dropping terms 0(e) with respect to terms retained in the shock layer 
it is seen that we recover the usual stream tube and normal momentum equations. 


Dividing (202) by the total pressure on the streamline and integrating, 
along the body surface from the stagnation point yields 


p(D/p t ~ 1 


( i ,j 

*1 + (2p (£)/Pu‘ 


which in the forebody as E, 0 approaches the stagnation limit 


LIM 

£ - 0 


p(5)/p t - i = o => 


p(€) 


-*• l 


In the afterbody as £ -*■ 00 the pressure approaches the vacuum limit (the stream tube 
"fault", unlimited decay of momenta) 


LIM 

5 ■+• o 


Elll Piil 


-»• 0 


r> 


So that an expected solution of the stream tube momentum equation subject 
to these limits near the stagnation point (variable area S^) would have the form 


5L.fi r o -> = - 2. 

p 5 C 5 » p 


(204) 
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or — = - Nd &N (£) 

,Pv . _ . . C N (20 

or (— J _«* 1 - Constant * £ 

v P'y = 0 
T 

Equation (204) is clearly inappropriate for the stream tube at large £ . 
However, as noted, an asymptotic limit pressure must be .applied to insure that 
the stream tube expansion limit is not exceeded. 

The dependence on curvature may be deduced from an integration of the y 
momentum equation. Equation (203) following division by the total pressure. 



dy 


6 


icpu 2 


0 pu z +2p 


dy 


(205) 


For £ *>■ °° pu -*• °° 


and 

T 


= f (y) 



(206) 


Hence p/p^ may be expected to approach a constant value along a given nor' 
mal. The strong entropy layer effect on pressure is absorbed by using the ratio 
of the local pressure to the "local" streamtube total pressure. 


For £ -»■ 0 on the other hand, pu 2 0 and from (205) 





(207) 


Or the total pressure variation in the forebody is linearly dependent on 
the curvature variation. In our simplification it is assumed to vary proportional 
to y. 


The present results correlate in the form implied by Equation (204) . Along 
the body streamline we find 


(-2-) _ = 1 - 0.247 (£) 0,328 (208) 

P T y = 0 
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For (J-) 2 (-£-) 


P T P T Limit 


(208) 


While 


tfi 


Ain z a + 0.11 4ina For 15° < a < 53.5° 


(209) 


T Limit 


In the shock layer (0 < y < 6) the presented results are adjusted to 
correlate in accordance with Equation (207) 


*T y > 0 *T y = 0 *- w o _ 


( 210 ) 


Here 


K . = K (maximum) = k (X = 0) 
w q w w 

The relative invariance of the streamtube pressure correlations to changes 
in reentry trajectory conditions is illustrated in the accompanying correlation 
plots Figures 7-10, 7-11 and 7-12 for the 7 cases (trajectory points) previously 
listed in Table 6-1. On these plots appear the correlation calculations developed 
by the procedures introduced here plotted against the mean values taken from the 
references previously listed. 

Tables 7-1, 7-2, and 7-3 summarize the numerics for the 7 angle of attack 
shock layer flow situations treated in the inviscid analysis. The streamtube 
pressure distributions, in terms of p/p^ shown in Table 7-2 are virtually in- 
variant over the trajectory range considered in this analysis however Mach number 
and angle of attack effects are present in the pressures, Pg/P T via the computed 
local conditions behind the shock wave. These initial conditions behind the shock 
wave are presented in Table 7 ” for each trajectory point. 


7-14 















7-17 


ACTUAL. VALUE. 

FkSiUtSe. *7-12 6Tt26\MT'J5e P(2eS5»U^S COieE£LLAT(CM 

^TeeAMLiue s am d s> 

















































































































































































































































k 


\ 


s 

M 


i 



W 


«o 

a. 

(Jl (J« O' ^ 

* *H i*> ® <0 <*> 'O S? 

<ji ® r* in **■ 

*■4000000© 

at o •-* * m o 

m •* o# « r- 

® ® ® r* \o *» 

*4 0 0 0 0 0 0 


S 

H 

— 

B 

gM 

ft 

M 

a 

p> h h n O' n « 

m ^ h if r* o' h 

° ^ ^ r< r-* m V 

rH *-« m V 

CD 

• 

O' 

m M Q o 1 * 

m * o o * ** 



Q PI CD PI O' O 

H N « 

a\ 

• 

3 

s 

§ 0 
r 

3 

sf^ouo^onu'i 

8dinr*oatOinin 
u «*♦ H ot <*» V 

S 

1 

Streamline 

7.0 

9.0 
11.0 
15.2 
20.0 
35.0 
46.5 

1 






«T 

a 

r~ 

• m «n ® O' in *-h 

<»l 01 O 0 U1 TJ r-» cM 

® O' O' TO r- >£> i/i in 

(Hooocoo do 

1. 

0.929 

0.910 

0.841 

0.764 

j.627 

0.514 

0.514 

1 

Q rH <N (A »*1 P^ 

® O « r- m m 

® 0 » ® ® m *n 

<-i o o o o o o 

1 

8 

M 

a 

iO 

m ® r- O o o OO 

pi hi «p « « v o -* 



OrHn # ® «r o *6 

•i n « 
at 

m 

in r-* in ® o o o 

w 'O m m « o o 

• 

r4 M IT ® ^ O ® 

rt pi a 
lO 

m 

<j\ r* \0 O c- O 

»D f O Cl H 

o * • • • • • 

«■* r> in <e O at 

1-4 <*> *T 

O' 

4 

M 

3 

§ 0 
S** 

8* 

P 

a 

|«»oonoOMc»a» 

# o H PI r! n r- m N ® 

u r-( m V 

4» 

u 

C 

5 

^OOMOOfJ^ut 

UHlNPli/lH^NlD 

M w* <n v 

** 
in 

0 

s 

^or-iOOM^at 

o ni r> k i r* m >1 -0 

3 * "» * 

(A 



7 


] 9 
















TABLE 7< 


1 


I 

I 








7 - 23 




i 



























i 


7.4 CHEMICAL RELAXATION ALONG STREAMLINES 

Using the pressure distributions given in Table 7-2 , the initial conditions 
given in Table 7-3 and an initially undissociated gas; the streamtube code was 
used to determine the thermochemical state along the streamlines. The chemical re- 
action data given in Column 1 of Table 5.7 was used since the comparisons of Sec- 
tion 5 indicate this set to be the more reactive of the two examined, thus pre- 
dicted heating rates should be conservative ly high. Calculations were made for 
t = 250, 400,, 800, and 1400 seconds and typical results for t = 800 seconds are 
shown in Figures 7-13 through 7-17. The mass flow represented by each of the stream- 
lines, was calculated from Table 7-1 and compared with the boundary layer flow rates 
shown in Figure 7-4. The 'match points' are shown in Figures 7-13 through 7-14 and 
therefore specify the boundary layer edge conditions at various locations on the 
body. Blending of these solutions with the viscous shock layer r ults (Section 
7.2) was achieved by plotting shock layer entropy (at constant n) as a function 
of X and selecting the entropy which blended best with the predicted downstream 
boundary layer edge entropies. This is shown in Figure 7-13; then other variables 
enthalpy and mass fraction, are taken at this same value of ri and compared with 
predicted downstream values. As shown in Figures 7-14 through 7-17, the tran- 
sition .rom shock layer values to predicted downstream values is surprisingly good 
except in the case of NO mass fraction where due to an overshoot behavior, small 
changes in "matchpoint" location x can result in large changes in mass fraction. 
Similar results were obtained for the t - 400 and 1400 second cases; however a 
significant interpretation difficulty :as noted in the t = 250 second case which 
is at a high altitude and therefore a low density. At these entry conditions, the 
shock layer is fully viscous so that the validity of a boundary layer analysis is 
questionable although shock layer solutions would still be valid. This question 
will be deferred to Section 7.6. 

From the above procedure, edge values of entropy, enthalpy and specie mass 
fractions for t = 400, 800 and 1400 were obtained and are shown in Figures 7-18 
through 7-22 as functions of the boundary layer coordinate, S. These were the 
distributions used in subsequent nonequilibrium boundary analyses. 
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7.5 


HEAT TRANSFER ON SHUTTLE VEHICLE 


Using the edge conditions given in Section 7.4 and the surface catalycities 
given in Section 5.6, nonequilibrium heat transfer rates to the pitch plane of the 
RI shuttle vehicle were calculated. The results are shown in Figure 7-23 along with 
the values predicted with equilibrium BLIMP and an isentropic edge expansion. 

Only the 400, 800 and 1200 second cases are shown since the question of the vali- 
dity of boundary layer assumptions at the altitude corresponding to t = 250 
seconds has not yet been resolved. The nonequilibrium assumptions are seen to 
reach a peak at a station slightly removed from the stagnation point whereas the 
equilibrium, isentropic edge solutions place peak heating at the stagnation point. 
This shift is due to entropy layer effects rather that chemistry effects as shown 
in Figure 7-24 where solutions using various assumptions are shown. 

From Figure 7-23 it is seen that some small benefit is attained by the 
reduced catalycity at the stagnation point however there is an apparent penalty 
downstream. Again the discrepancy is due to entropy layer effects as shown in 
Figure 7-24 where nonequilibrium and equilibrium solutions with entropy layer are 
compared . 

It appears from Figure 7-24 that the benefits of the slightly reduced 
catalycities derived in Section 5.6 have only small benefits however fully non- 
catalytic walls can reduce heating rates by 25 -50% throughout the length of the 
vehicle. 

One other possible anomoly appears in Figure 7-24. The semicatalytic and 
equilibrium entropy layer curves intersect at about s = 45 ft. and is mainly a 
chemistry effect. In the nonequilibrium case the boundary layer was close to fro- 
zen so that equilibrium assumptions would result in thicker boundary layers since 
homogeneous recombination behaves like a source. Thicker boundary layers result 
in reduced heat transfer which beyond s ~ 45 ft. is apparently greater than the 
effect of reduced surface recombination for the assumed surface catalycities. 

The accuracy of the results shown in Figure 7-23 can be improved by iter- 
ation on the edge boundary conditions. The existing solution can be used as a 
better approximation to the boundary layer mass flow for matching to the invscid 
solution. Note that, unless displacement effects are large (which would be the 
case at very high altitudes) , the inviscid solutions do not change. 
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7.6 


APPLICABILITY OF BOUNDARY LAYER ASSUMPTIONS AT LOW DENSITIES 


The stagnation pressure for the t = 250 seconds case is an order of magni- 
tude less than that at t = 400 seconds, namely 0.0014 compared to 0.018 atm. At 
t * 400 seconds a matching of shock layer and boundary layer solutions similar to 
that shown in Figure 7-13 through 7-17 was obtained even though the viscous zone 
occupied a significant portion of the shock layer. However, at t = 250 seconds, 
boundary layer predictions with BLIMP and BLIMP/KINET show a thickness greater 
than the predicted shock layer thickness. A comparison of shock layer and 
boundary layer solutions are shown in Figures 7-25, 7-26 and 7-27. Shock layer 
assumptions predict a stand-off distance of about 0.16 feet whereas various 
boundary conditions imposed on the boundary layer (see Figure 7-25) resulted in 
predicted boundary layer thicknesses between 0.2 and 0.45. Since this is not 
physically possible one must conclude that either the shock layer predictions 
or the boundary layer predictions are invalid. The latter is most likely since 
the predicted shock stand-off distances are in general agreement with the pre- 
dictions of Reference 76 f r a hypersonic density ratio of 1/10. Further, as 
shown in the shock layer curve in Figure 7-25, the influence of the wall extends 
throughout the shock layer; that is, there is no inviscid region. From 
Figure 7-27 it is apparent that the shock layer results cannot be approximated 
by frozen boundary layer edge conditions (curves A and D) . 

Although near the stagnation point, boundary layer assumptions are not 
valid, thin shock layer assumptions can lead to a better, though rigorously 
not completely valid, predictions. However downstream from this region the 
shock layer is not thin so that these predictions become invalid and still further 
downstream, conditions may be such that boundary layer assumptions are once 
again valid. Thus, although solutions can be obtained with boundary layer 
assumptions, one would be hesitant to place much confidence in them since initial 
conditions and boundary layer edge conditions cannot be adequately defined. 
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SECTION 8 


SUMMARY AND CONCLUSIONS 

A nonequilibrium boundary layer code has been developed which retains all 
of the boundary condition generalities of the equilibrium BLIMP code. Both 
laminar and turbulent boundary layers are permitted however no attempts were 
made to resolve the question about the validity of standard chemical kinetics 
models in turbulent flows. The code BLIMP/KINFT was used for extensive 
studies on the boundary layer and shock layer behavior in the stagnation region 
of both shuttle and RV size vehicles. 

A procedure has been developed to approximate nonequilibrium and non- 
isentropic the rmc chemical states at the edge of a boundary layer on the pitch 
plane of typical shuttle vehicles. 

The catalycity of typical shuttle TPS materials in dissociated air 
was estimated from a sampling of arc jet data. This data and the boundary layer 
edge conditions were used with the BLIMP/KINET code to calculate shuttle heating 
rates in laminar and turbulent flows. Some conclusions from this investigation 
are given below 

1) The catalycity of typical TPS materials such as LI -900, LTV coated 
carbon-carbon and coated columbium is estimated to be about 1020 cm/ 
second for 0 recombination and about 312 cm/second for N recombina- 
tion at enthalpies between 2000 and 1.4000 Btu/lbm and pressure of about 
0.01 atmospheres. 

2) Entropy layer effects are essential whether equilibrium or nonequilib- 
rium effects are considered. 

3) Noncatalytic walls can reduce heating rates on the shuttle vehicle by 
25 - 50% however, for the catalycities estimated from arc jat data only 
minimal nonequilibrium effects are predicted. 
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4) For the portion of the entry trajectory betw*>r. t = 250 ax'd t ** 1400 
seconds of pressure ratio, the distribution (P/l?^) along streamlines 
is not sensitive to altitude, velocity and small angle of attack 
changes. These effects are implicit in the initial conditions behind 
the shock wave. 

5) In most cases nonequilibrium chemistry calculations do not require any 
more nodal points t.*an the comparable equilibrium calculations . 
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APPENDIX A 


INPUT DATA FOR SURFACE CATALYZED REACTIONS 


Surface reaction rates are calculated as if heterogeneous chemistry were 
of the same form as homogeneous chemistry, i.e., for the recombination of oxygen 
atoms, the net surface reaction is 


20 -° 2 


(A-l) 


The molar production rate is then 


(MJ = Mp* -f- p 0 ) 

P 2 


(A-2) 


The term Po 2 /K will be small for nonequilibrium conditions so that 


(M q ) - k £ ( P J) 


(A-3) 


This relationship indicates that the catalysis reaction is second order due to 
the form of Equation (A-l) . However there is some evidence that the reaction 
is first order or possibly something intermediate since catalysis is not really 
a one step reaction. Any reaction order can be approximated by rewriting Equa- 
tion (A-l), For instance, to specify surface recombination as a first order reac- 
tion, one can write 


0 


I°2 


(A— 4) 


Then Equation (A-3) becomes 


«o> = W 


lb-moles 

ft 2 sec 


(A-5) 



In order to express k^ in terms of surface catalycity, one can write 


p k (SP.) 
w w i w 



m 


w 


^ (M i } 


lb-moles 

ft 2 sec 


) 


(A-6) 


For perfect gases 


P (SP.) = 

w i w 


pA 

RT 


(A-7) 


so that 


p (SP. ) RT 
' W 1 w 

k f th ± 


p k (SP.) 

W W 1 w 



lb-moles 
ft 2 sec-atm 

Thus for the first order reaction (A-4) k and k- are related by 

* W X 

k f = (0.025 k w ) i 

where T (°K) and k w (cm/sec) . In terms of input variable 

FKF = 0.025 k 

w 

POW = -1.0 
EAK = 0.0 



(A-8) 


Similar relationships could be obtained for other reaction orders. Note that in 
the above formulation a dependence of k w on temperature could be accounted for 
in the form of 


A-2 



k 

w 



„ T n e (-E/RT) 


where k w , is some reference value. Then 
w ref 


PKF = 0.025 k 

w * 
ref 


POW = -1.0 + n 
EAK = E 



